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Abstract 


In this thesis we consider the solution to dynamical ice flow equations using a com- 
bination of a moving mesh method and data assimilation. We show that by moving 
the mesh the approximation to the ice thickness profile is improved, and the location 
of the domain boundary is significantly better estimated. The method used is derived 
by utilising a relative mass conservation principle to define a net deformation velocity 
comprising of the internal diffusion of ice and the effect of accumulation or ablation. 
We use a finite difference numerical approximation in one-dimension and a finite ele- 
ment approximation in two-dimensions to demonstrate the ability of the methods to 
simulate different aspects of ice flow. In particular we focus on the accurate represen- 
tion of the moving front of the glacier without the need for an interpolation procedure. 
Results are shown to compare favourably to exact, steady state solutions, while demon- 
strating improvements over traditional fixed grid methods. 

The impact of the internal diffusion of the ice on the movement of the glacier front is 
analysed, and a condition on the local profile near the boundary is constructed to de- 
termine when the front is moving as a result of diffusion rather than the accumulation 
or ablation. 

We utilise the technique of data assimilation to combine the moving mesh method with 
observational information to get a statistically best estimate of the ice thickness profile. 
In a moving mesh environment there are differences to the scheme that we detail, in 
both one and two dimensions. 

We introduce an extension to our data assimilation scheme to directly include the nu- 
merical mesh within the update. This allows for the potential inclusion of observations 
of key features such as the location of the boundary. We demonstrate the improvement 
that this extension has on our prediction of the domain in one dimension and discuss 


the challenges encountered when applying this extension to two dimensions. 
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Chapter 1 


Introduction 


Advances in technology and the ability to reach the far ends of the Earth has allowed 
the study of the polar regions to accelerate in recent times. Within the polar regions, 
as well as alpine areas, there exist large bodies of ice that for the most part appear sta- 
tionary. In fact the ice diffuses slowly under its own weight, and its movement impacts 
upon both the natural environment and human civilisation. 

The state of ice on planet Earth has been the focus of much attention in recent 
times, particularly the rate at which the ice mass is decreasing. As with many natural 
features the impact of ice on the global climate is greater than in just the polar or 
alpine regions in which it exists. Ice stores around 80% of the fresh water found on the 
planet, which were it all to melt, would greatly impact the salinity of the ocean water 
as well as cause the much publicised sea level rise (~ 65m if all ice melted). This would 
impact ocean circulations such as the Gulf stream which could ultimately bring long 
term changes to climate systems [75]. In addition the albedo of ice is high, between 0.5 
and 0.9, which acts to reflect large amounts of incoming solar radiation [42]. On a more 
local scale a loss in ice mass can cause natural disasters such as floods and avalanches. 

Moreover, during the various ice ages, glaciers and ice sheets expanded to claim 
large areas of the planet resulting in loss of life and, for some species, extinction. In 
the past the advance of ice has been likened to life threatening monsters such as the 


ice dragon drawn by HG Willinck (Figure 1.1). The glacier pictured (taken from [13]), 
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the Mer de Glace is no longer visible from the same location due to an extended period 


of retreat. 


Figure 1.1: The Mer de Glace in the French Alps depicted as an ice dragon [13] 


The recent news headlines claiming that all ice will melt within the next couple of 
generations are generally considered unlikely; however, the most recent assessment re- 
port of the IPCC (Intergovernmental Panel on Climate Change) in 2007 suggested that 
mathematical models are not yet sufficient to make credible predictions [96]. With this 
in mind much work has been devoted to increasing the accuracy of these mathematical 
models to reduce their uncertainity and provide better predictions for the future. These 
uncertainties in the models arise from two primary sources; 

Firstly the dynamical response of ice to varying stress levels is extremely complex. 
While ice is generally considered to be viscous, it is highly non-linear, occasionally 
anisotropic and on some short time scales viscoelastic [51]. In addition there is the 
presence of debris, water and air mixed in between the ice. 


Secondly the amount of data available is severely limited. Surface and satellite ob- 
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servations can give details about the upper layers of ice, but the state of deeper ice 
is harder to determine and is crucial to the flow of the entire body. For example the 
presence of meltwater near the base can cause basal sliding and thus speed up flow, 
which is currently an unresolved issue even if the location of any meltwater is known. 

It is only recently that dynamical ice flow models have started to be included within 
global climate models. Part of the reason is the large uncertainties contained within 
the models, whilst also being computationally expensive, requiring a higher resolution 
than is present in current climate models. One potential solution is to use a simpler 
dynamical model which requires less computational time, though this invariably in- 
creases the uncertainty. Alternatively, ice sheet and climate models could be operated 
separately and coupled together. Generally this option would involve a coarse coupling, 
using ocean and atmospheric conditions from a climate model as forcing terms in the 
ice sheet model, and the ice cover extent and water outflow in the climate model [51]. 
Neither is a particularly elegant solution to an open problem. 

An alternative to increasing the resolution is to consider adaptive mesh techniques, 
which allow better approximations to be made in the relevant areas. Several current 
methods utilise mesh refinement to increase the resolution in certain areas which work 
well provided the domain is reasonably well understood and the areas that require 
higher resolution are approximately known and unchanging. There are also a number 
of methods which move the computational mesh nodes to the required areas and have 
been shown to be more robust than non-adaptive methods [105]. Currently these move 
the nodes independently of time and are controlled by a monitor function which can be 
difficult to define [22]. In some other fields the movement of the mesh nodes is defined 
in terms of a time dependent velocity, which provides a clearer definition of the monitor 
function and allows the mesh to be influenced by its previous position(s). Currently 
this approach has not been applied to the dynamical ice flow equations. 

The amount of information from satellites, aircraft and ground based measurements 
of ice sheets and glaciers is increasing. The glaciology community is beginning to make 
use of these observations in conjunction with numerical models using the techniques 


of data assimilation to produce the best possible representation. The atmospheric and 
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oceanic communities have been developing these techniques for many years, but the 
application to glaciology is in its infancy. In addition the application of data assimi- 
lation has largely been restricted to fixed grid methods, naturally since these are the 
most prominent methods used in operation. However it has been shown that adaptive 
methods are useful tools in the right situation [99], and at some point there will be a 
requirement to combine numerical adaptive methods with data assimilation. 
Therefore there appears to be a gap in both the use of moving meshes and data 
assimilation applied to the field of glaciology. There also has been little development 
of data assimilation methods on moving meshes regardless of application. It is this gap 
that we seek to begin filling. We now lay out our aims for this thesis, followed by an 


overview of the key results and a detailed chapter-by-chapter plan. 


1.1 Aims 


In this thesis our main objectives are to: 


e Develop a moving mesh method to equations governing dynamical ice flow. 


e Use data assimilation to combine observations with a moving mesh model of 


dynamical ice flow. 
In particular we aim: 


1. To apply the conservation of mass fraction moving mesh method of Baines, Hub- 
bard and Jimack [8], henceforth known as BHJ, to the one-dimensional shallow 


ice approximation equation used within glaciology. 
2. To analyse the impact of internal glacier flow on the movement of the ice boundary. 


3. To use a sequential 3D-Var data assimilation scheme in conjuction with the mov- 


ing mesh to improve the predicted ice profile. 
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4. To build an extended scheme capable of including observations of positional fea- 
tures within the data assimilation method to improve the representation of the 


domain. 


5. To extend the moving mesh method and the 3D-Var scheme to two horizontal 
dimensions and demonstrate the difficulties encountered when considering this 


more complex scenario. 


1.2 Key Results Achieved 


In this thesis we show that: 


e Moving mesh methods are an effective tool to model the shallow ice equations 


quickly and accurately. 
e Data assimilation can be performed on a moving mesh domain. 


e Including the numerical mesh within the data assimilation scheme can improve 


the accuracy of the domain prediction in 1D. 
Specifically we show that: 


1. Using a moving mesh scheme based upon conserving fractional mass to simulate 
a 1D shallow ice flowline model can accurately simulate the key features of the 


ice sheet with minimal computational cost. 


2. The impact of the flow velocity at the boundary is dependent on the local profile 
of the ice thickness, with an infinitely steep gradient required for the boundary 


to flow. 


3. Results compare favourably to both the exact solution and fixed grid approxima- 


tions to the European Ice Sheet Modelling INiTiative (EISMINT) experiment. 


4. A sequential 3D-Var data assimilation scheme can incorporate observations into 


a moving mesh dynamical model. 
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5. With the chosen moving mesh method a recalculation of the total mass is required 


when assimilating to ensure mass conservation. 


6. By extending the data assimilation scheme in 1D to include the numerical mesh 


we can find an improved solution of the mesh, enabling greater accuracy. 
7. The moving mesh scheme readily extends to two horizontal dimensions. 


8. The 3D-Var scheme can be applied to a 2D moving mesh to improve the numerical 


approximation. 


1.3 Thesis Outline 


In chapter 2 we begin by introducing the field of glaciology. We highlight the main 
driving forces behind the flow of ice and introduce the full Stokes model for simulating 
this flow. We then discuss how this is simplified to the shallow ice approximation model 


relevant to the work in this thesis. 


We then give an introduction to adaptive mesh modelling in chapter 3, describ- 
ing different ways a mesh may be adapted to improve the numerical approximation to 
the solution. In particular we focus on a moving mesh scheme, the relative mass con- 
servation method of BHJ [8] that we will be applying to the shallow ice model in this 


thesis. We then discuss the current use of adaptive meshes within the field of glaciology. 


Chapter 4 begins with a general introduction to data assimilation and some of the 
algorithms used operationally. We discuss the reasons for choosing the 3D-Var sequen- 
tial method in this thesis over the other methods. We then look at the use of data 
assimilation within glaciology and the types of data that are available. Finally we dis- 


cuss the attempts at using data assimilation on adaptive meshes. 


In chapter 5 we analyse the behaviour of the internal flow velocity and its effect 
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on the boundary, providing conditions on the local glacier profile required to induce 
movement at the glacier front. We then apply the BHJ scheme to the shallow ice 
approximation. The scheme is discrestised to form a numerical moving mesh method 
using finite differences as a means of approximation. We use a test scenario to illustrate 
the ability of the method to replicate aspects of glacier flow, before using the radially 
symmetric EISMINT problem to assess the benefits of the approach over existing meth- 


ods. 


We apply the 3D-Var data assimilation algorithm to the moving mesh method in 
chapter 6. We determine the required components of the algorithm and the modifica- 
tions needed to apply it to a moving mesh model instead of a fixed grid. This is tested 
using a twin experiment applied to the test scenario introduced in the previous chapter. 
We propose an extension to the algorithm to include the numerical mesh within the 
data assimilation algorithm, allowing the mesh do be updated at the assimilation time 
and the potential to include observations of positional features. The changes required 
to the components of the algorithm are examined, before we demonstrate the impact 


on the solution. 


In chapter 7 we extend the moving mesh method to two horizontal dimensions. We 
rewrite the method in a weak formulation, so that we may discretise the method using 
a finite element approximation. Again we use a test scenario to assess the movement 


of the mesh, before testing against the EISMINT scenario. 


We extend the 3D-Var algorithm to be applied to the 2D moving mesh method in 
chapter 8. We determine the components of the algorithm for use with a 2D moving 
finite element model. Using a twin experiment we test the 3D-Var scheme, demonstrat- 
ing an improvemed estimation to the true ice thickness profile. We then discuss the 
extended algorithm and the options available to apply the algorithm, with the benefits 
and difficulties of each. 
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Finally, in chapter 9 we summarise the conclusions made throughout this thesis and 
discuss the limitations of the methods introduced. We highlight potential extensions 


for future development of these methods. 


Chapter 2 


Glaciology 


In this chapter we examine the rheology of ice and the way this has been implemented 
into a definition for the interior deformation for use within mathematical models. We 
examine the full Stokes equations in the context of dynamical ice flow and the sim- 
plification to the Shallow Ice Approximation (SIA). We detail the assumptions that 
are made in this assumption order to compute forecasts, with a detailed look at the 


equations that will be used in this thesis. 


2.1 Introduction 


Glaciers and Ice Sheets are subjects of much speculation over their contribution towards 
sea level rise as ice melts in response to climate change. As ice is removed from the 
glacier system the behaviour of the remaining ice can change dynamically as a response 
to this melting. It is therefore important to understand the processes that drive dy- 
namical glacier flow in order to accurately simulate changes that may occur. 

At this stage it is worth introducing a few definitions to describe the different fea- 


tures of grounded ice: 


e Ice Sheet - A glacier of considerable thickness and more than 50,000 sq. km 


in area, generally slow moving across a large scale. However, they tend to feed 
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many outlet glaciers. Currently only Antarctica and Greenland are classed as ice 


sheets. 


e Ice Cap - Ice caps are a dome-shaped cover of perennial snow and ice. They are 


classified as small ice sheets, generally less than 50,000 sq. km in size. 


e Ice Stream - Ice streams are channelised glaciers that flow more rapidly than 


the surrounding body of ice. 


e Alpine Glacier - A glacier occupying a valley, usually high in mountainous 


terrain. 


e Outlet Glacier - Glaciers which emerge from larger ice sheets but are constrained 


like alpine glaciers. 
e Glacier Front - The moving boundary of ice that advances or retreats. 


e Ice Divide - A topographic separation of either multiple glaciers or a glacier /rock 


interface where there is zero horizontal velocity. 


There exists a variety of mathematical models to describe the movement of ice, ranging 
in complexity. We now present the key driving forces behind the movement of ice 
before presenting a selection of these models, concentrating solely on ones which apply 


to grounded ice. 


2.2 Ice Rheology 


For the most part ice sheets and glaciers are made up from Polycrystalline ice, where 
the composing crystals vary in orientation and size. This is due to the formation of 
ice from snowflakes landing with randomly distributed grains, leading to the overall 
behaviour of ice being treated as isotropic [48]. The behaviour of ice is not always 
isotropic and work has been done on implementing anisotropic behaviour into models 


developed for isotropic ice (see e.g. [73]). In this work we shall retain isotropic ice 
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behaviour. 

Under low levels of stress ice behaves in a viscoelastic manner, whereas under larger 
levels of stress the elastic effect is negligible and the ice behaves like a non-linear, non- 
newtonian fluid. In general ice is treated as viscous with the elastic behaviour ignored. 

In an ideal world a mathematical equation linking together the key components of 
stress, rate of strain and ice properties such as temperature would exist to describe dy- 
namical ice flow. However this has not yet been discovered, if indeed such a relationship 
exists at all. In order to gain some understanding of how ice moves, some simplifica- 
tions need to be made. If it is assumed that ice is incompressible, a not unreasonable 
assumption is most cases, then the strain rate (€) is proportional to deviatoric stress 
(o’) [82]. This provides a common constitutive relationship acting in the cartesian 
directions x, y, 2: 


Eig OC O35 eae nee (2.1) 


The most commonly used relationship is Nye’s generalisation of Glen’s flow law 


(40; 83] 
é 1-1/n 
5 e / 


where A and n are flow law parameters and €, is the effective strain rate. Most dynami- 
cal modelling of ice flow is based upon extracting the velocity field from this relationship, 
which works well but is by no means a complete description of ice flow. 

The two flow law parameters are unknown quantities which have empirical estimates 
from laboratory experiments. Generally it is accepted that the flow law exponent n = 3 
is sufficient in the majority of cases [108; 44]. For the rate factor A, a modification is 


made to the Arrhenius equation to fit the laboratory and field data [44]: 


Q 3C 


art mF: (2.3) 


A= Agexp |— 


where T' is the temperature of the ice and the remaining variables are constant with 
values are given in Table 2.1. The flow law parameter in Eq. (2.3) is dependent on the 
temperature of the ice and has a log-linear relationship, except close to the melting 


point. The variability in A is quite drastic and has a large influence on the speed of 
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Symbol Value Definition 

Ao 9.302.10-2Pa-7yr-! Normalising factor 

Q 7.88.J(mol)~! Activation energy 

R 8.314] (molk)~! Gas constant 

To 273.39K Reference temperature 

C 0.49836.K* Empirical constant 

k; i rg Unitless constant 


Table 2.1: Table of typical values for Hookes modified Arrhenius equation. [13] 


Rate factor A (kPa? yt) 


-50 -25 0 


Temperature (°C) 


Figure 2.1: The rate factor A, determined by [27; 44]. The value of A, and as a result the 


speed at which ice diffuses varies dramatically as the temperature of the ice changes. (Image 


from [102]) 


the ice flow (see Figure 2.1). 


Glen’s flow law is not the only relationship between stress and strain. An alternative 
is the Smith-Morland flow law [94], which solves the problem of infinite viscosity when 
o, — 0 that can occur within Glen’s flow law. According to [94] it also fits the laboratory 
experimental data better then Glen’s flow law, though since is not widely used we shall 


retain the use of Glen’s flow law in this thesis. There are also flow laws for dealing with 
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anisotropic ice (see [5; 38]), which we will not consider here. 


2.3 Mass Continuity 


The mass continuity of ice may be derived by considering a column of ice stretching 


from the base to the surface (see Figure 2.2). 


Fant 


= F base 


Figure 2.2: Mass changes in a column of ice (Modified from [10]) 


The flux of ice entering this column from the rest of the glacier is given by: 


F;,, = inflow in x-direction + inflow in y-direction (2.4) 


= hu* Ay + hu’ Az. (2.5) 


where h = h(x,t) is the thickness (or height) of the ice and (u”,u¥) is the ice velocity 


in the Cartesian directions. The ice leaving the column is expressed by: 


Fout = outflow in x-direction + outflow in y-direction (2.6) 


= (h + Ah)(u® + Aut) Ay + (h + Ah) (uY + Au?) Ac. (2.7) 
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Mass can also differ in the column due to accumulation/ablation at the surface and 


melting/refreezing at the base, 


Eig pig ey (2.8) 


Pigss = ™mpArcAy, (2.9) 


where m, and my, is the mass flux at the surface and base respectively and can be 


positive or negative. The change in volume therefore is: 


OV Oh 
cBE. = ap Ot AY = E; =e Lauer ooo bane (2.10) 
= —(u°Ah + hAu*)Ay — (w¥Ah + hAu”) Ar 


+m,Ardy + mAzAy + O(A?) (2.11) 


Ignoring higher order terms, we divide through by AzAy, leaving the change in ice 
thickness as Ax, Ay — 0: 


Oh wAh+hAu® wu¥Ah+ hAu 
a Ae - i +m, + mp (2.12) 
x y 
_ a )_ — cee (2.13) 
= —V.(hu) +m (2.14) 


Eq. (2.14) is an equation of mass conservation as defined in [72], incorporating 


diffusivity and an external source term. 


2.4 Models 


The mass continuity equation (Eq. (2.14)) can describe the movement of ice thickness 
over the glacial domain, provided we are given the source term, m, and a description of 
the flow velocity, u. The source term is provided from external information, while the 


internal flow velocity requires further information. The most complete mathematical 
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description of dynamical ice flow comes from the set of full Stokes equations, which 
we will detail here. Due to the complexity of these equations many approximations to 
these equations have been made under different assumptions. We will also detail here 
the most commonly used simplification, the Shallow Ice Approximation (SIA) which 


we use throughout this thesis. 


2.4.1 Full Stokes model 


The full Stokes equations for ice sheets begin with the equations of motion [48]: 


O 
a +V.(pu) = 0, (2.15) 
Ou , 
P\ > +(u.V)u } =—Vp+ V.o' + pg, (2.16) 


where p is density, u is velocity, p is pressure, o’ the deviatoric stress tensor and g body 
forces. Ice is a slow moving body, which means the acceleration terms in the momentum 
equation are negligble. In addition there is very little change of density with depth or 
temperature, as the layers of snow and firn (old snow) are shallow compared to the 
main body of ice; in other words the body of ice may be treated as incompressible. 
Lastly the only body force acting on ice is gravity. Therefore the equations of motion 


Eqs. (2.15) and (2.16) can be written in the form 
V.u=0, (2.17) 
Vp =V.o' + pg. (2.18) 
Development of full Stokes models is still in its infancy. Currently a typical Stokes 
model will take a long time to run on a serial computer [56; 55]; however, there is 


an open source code available which is designed to run on parallel systems to reduce 


computational time [1]. 


2.4.2 Shallow Ice Approximation (SIA) 


One of the most widely used simplifications of the Stokes equations is the Shallow Ice 


Approximation (SIA). It is a zero-order model and while it is generally considered to 
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be a coarse approximation; computationally it is very cheap. The SIA is based upon 
three primary assumptions [48]: 

Firstly the ice sheet is considerably wider than it is thick, and the gradients are small. 
Le. 


Az << Ag, Ay: (2.19) 


This is a reasonable assumption, as Az, Ay ~ O(100km+) and Az ~ O(1km). Sec- 


ondly it is assumed that the flow of ice is primarily in the horizontal directions, with 
I (2.20) 


Again this is a reasonable assumption for many scenarios; for example the Vostok ice 
core showed it took ~ 400,000 years for ice to travel ~ 3km in depth [87]. This means 
that the vertical strain rate €,, is set to zero. Thirdly, the assumption is that the change 
in stress and consequently, strain, is negligible in the horizontal directions compared to 


vertically. i.e. 
OSes. . Otiy 


aye: (2.21) 
Applying these three assumptions leaves the following equations 
= = (= iz =) . (2.22) 
oP SOs (2.23) 
" = ou (2.24) 
*p = pg. (2.25) 


Since density is constant, we may depth average Eqs. (2.22) to (2.25) to further reduce 
the number of variables to find the flow velocity u, using Glens flow law (Eq. (2.2)) to 
give 

2A 


So ph gh pete il WS 2.26 
a ong |Vs| 8 (2.26) 


where s is the surface elevation. The combination of Eqs. (2.14) and (2.26) provides a 
full description of the shallow ice approximation [48; 77]. Further details on how to get 


from the balance equations to Eq. (2.26) can be found in [48; 102; 13]. 
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2.4.3 Other Models 


There are many models which sit within the spectrum between full Stokes and the SIA. 
Most of these sit within a class of models known as Longitudinal Stress Approximations. 
The general aim here is to start with the SIA and re-introduce horizontal gradients 
of longitudinal stress into the momentum equations, which result in changes to the 
constitutive relationship Eq. (2.2). The variations in the methods of including this 
stress and the approximations made provide a variety of models [43], for example one 
of the main models in this class computes the longitudinal stresses on various layers 
throughout the ice sheet [85; 14], as opposed to at the surface like most of the models 
in the class, e.g. [71]. While each of these provides a better approximation that the 
SIA, the increased level of complexity means there are fewer examples out there with 


which to test our numerical method. As such the SIA will be used throught this thesis. 


2.5 Summary 


In this chapter we examined the rheology of ice and the definitions of stress and strain. 
We showed how these concepts relate to form a flow law describing the deformation of 
ice under its own weight. 

We then derived the mass continuity equation by considering the fluxes into a column 
of ice to obtain the main equation that will be used in this thesis. 

Lastly, we detailed the full Stokes equations for conservation of momentum in glacier 
flow, before making a number of assumptions to simplifly these equations to the shallow 


ice approximation. 


Chapter 3 


Numerical Modelling 


In order to solve differential equations such as the shallow ice equation in glaciology we 
turn to numerical methods to approximate the solution of the mathematical equations. 
In this chapter we look at the different types of adaptive mesh modelling and the 
benefits they have over traditional fixed grid methods. We introduce a type of adaptive 
mesh model where the nodes of the mesh are moved in accordance to a prescribed 
requirement — that the relative mass is conserved throughout, to provide a dynamical 
evolution of the mesh. We then discuss the types of adaptivity that have been used 


within numerical models of ice sheets and define some terminology for use in this thesis. 


3.1 Introduction 


Many physical problems are evolutionary in nature and can be represented in mathemat- 
ical terms by time dependent partial differential equations (PDEs). Of these physical 
problems many are non-linear in nature, exhibit unpredictable behaviour and contain 
localised transient features which are difficult to approximate numerically. Tradition- 
ally, standard numerical solutions to these PDEs are achieved by taking the physical 
domain and solving on a computational one by dividing the domain into a set of grid 
points. The solution is then approximated on each grid point and pieced together to 


gain an overall picture across the domain. Evolutions of the approximate solution to 


18 


CHAPTER 3. NUMERICAL MODELLING 19 


the differential equations upon this domain then formulate solutions over space and 
time. 

The basic form of a numerical scheme is a so-called fixed grid method, which has a 
rigid regular structure that is generally well understood and straightforward to achieve. 
There are however limitations in this approach. For instance, in some areas of the 
domain there may be features that exhibit sharp changes in some quantity over small 
spatial scales, the most obvious being a steep physical gradient. The ability to rep- 
resent these features is limited by the distance between grid points. Since the grid is 
generally uniform, reducing the distance between grid points will decrease the spacing 
globally, which increases the resolution of the numerical solution and achieves greater 
accuracy. However increasing the resolution across the whole domain also increases the 
computational time, sometimes to an impractical level. 

Another limitation with fixed grid structures lies with problems where the domain 
in question alters in size. To resolve this the computational domain needs to be large 
enough to cover any change in size that potentially could occur. This leads to a num- 
ber of grid points effectively going unused, which uses more computational time. In 
addition the location of some features, such as the boundary, will require interpolation 
or extrapolation to resolve their location on the domain. 

With these limitations in mind a number of alternative methods exist that modify the 


structure of the numerical domain to reduce errors in the solution. 


3.2 Adaptive Meshes 


As an alternative to formulating a numerical solution on a rigid structured domain, 
it can be beneficial to adapt the grid and the numerical solution to suit the problem. 
These methods have existed for many years with the aim of creating irregular grids on 
which to numerically approximate a more accurate solution than a fixed grid method. 
Adaptive meshes provide an efficient simulation for many physical problems. Broadly 


speaking these adaptive schemes can be split into three main categories; 
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e h-refinement - The spatial mesh is refined in regions where higher accuracy, and 
thus increased resolution is required. This increases the number of nodes by 


subdividing the relevant mesh cells. 


e p-refinement - The order of numerical approximation is increased locally in areas 


where more accuracy is needed. 


e 7-refinement - Also known as moving mesh methods. The number of nodes re- 
mains constant but the location of these nodes are distributed across the grid to 


areas where they can provide the most information. 


In addition these types of refinement can also be combined to improve the robustness 
of the representations (see e.g. [6]). 

As a compromise between fixed grid and adaptive mesh methods, the h- and p- 
refinement techniques can be used at an initial time to provide greater accuracy in re- 
quired areas. These methods are still fixed grid methods with a predetermined variable 
resolution. They do however require prior knowledge, exact or estimated, in order to 
determine the location for the refinement, which will in general change with time. 

Instead it is preferable for adaptive meshes to respond to the evolution as part of the 
adaptive process such that the refinement is performed automatically. To implement 
such schemes a solution indicator is required to determine the adaption. There are 
many methods for defining the solution indicator: the most common is the use of error 
estimates; however, a variety of other techniques exist. 

The dynamical ice flow model presented in this thesis contains features where high 
resolution is required, most critically close to the location of the moving glacier front. 
In steady state the location of these features would be known and a pre-refined or 
h-refinement technique would be applicable. However glaciers are rarely in a steady 
state and the boundary is constantly moving over time. For a domain that is frequently 
changing shape and size, r-refinement techniques are well suited to provide an efficient 
approximation as the moving boundaries are implicitly defined. 

While the use of r-refinement methods tend to be less popular than other adaptive 


techniques, they can be incredibly useful in numerical approximations. In particular 
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Figure 3.1: Example mesh mapping from a uniform computational domain (left) to an 


irregular physical domain (right). [17] 


the time-dependent expanding or contracting nature of many problems can be better 
approximated without the need for extra grid points. It can also be demonstrated that, 
in one-dimensional problems at least, r-refinement methods can resolve sharp transient 
features without the need for an excessive increase in mesh nodes [74]. 

There exists a number of techniques for generating the nodal movement in moving 
mesh methods, which can be classified into two subcategories; location and velocity 
based methods [22]. 

In location based methods the mesh is redistributed at each time step by directly 
redefining the positions of the nodes. This generally involves mapping from a uniform 
computational grid to the r-refined physical grid, like the 2D example in Figure 3.1. 
The mappings for these methods are often found by minimising some functional related 
to a mapping indicator, sometimes known as a monitor function. One such indicator 
is the equidistribution principle (e.g. [46; 12]) which chooses some property, such as 
the arc length, to be distributed evenly over the physical domain. Another commonly 
used indicator is the energy of a harmonic mapping [100]. The construction of these 
mappings is straightforward in one dimension, but can be more difficult in higher di- 
mensions [8]. In addition they need to be converted to velocities for incorporation into 
time-dependent PDEs such as the mass continuity of ice, Eq. (2.14). 


Velocity based techniques solve for the time derivative of the nodal positions, known 
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as the mesh velocity. The mesh velocity field is then integrated to provide the new mesh. 
As these methods are similar to Lagrangian type methods used in fluid dynamics they 
suffer from some of the same deficiencies, such as mesh tangling [18]. 

In the work of Miller and Miller [76] the mesh movement is approximated by min- 
imising the PDE residual in a finite element framework. A deformation map method has 
also been used by a number of authors [66; 20]. Cao et al.[21] have used the geometric 
conservation law (GCL) to form a functional to find the velocity through minimisation. 

A similar approach uses a local conservation principle, a modified version of the 
GCL, to directly find discrete velocities [7], which we will now describe in more detail 


such that it may be applied to the glaciology problem in this thesis. 


3.3. A Relative Mass Conservation Method 


One approach for defining the velocity at any time interval is the BHJ method |7; 8). 
This method is applicable to any mathematical model that is a balance between the 
rate of change of a volumetric quantity, the flow of said quantity across the volume 
boundary (known as flux) and any potential source or sink, where there is a change to 
the overall mass of the domain. As such the monitor function is chosen to control the 
relative mass of a set of subdomains which cover the entire domain and may overlap. 


Consider a non-mass conserving PDE of the form 
h, = F(h) +m, (ool) 


where t is time, h = h(x,t) is a positive dependent variable such as mass, F is a 
nonlinear spatial operator containing a combination of h and its spatial derivatives, and 
m is a source function independent of h. The total mass may be found by integrating 


the dependent variable over the entire domain and is defined by 


_ : hdx. (3.2) 
Q(t) 


Since the PDE, Eq. (3.1), is non-mass conserving the variable @ varies over time, so 


a relative mass conservation principle is introduced by normalising the mass in any 
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Figure 3.2: Relative mass conservation concept for a velocity based moving mesh method. 
The mass of w, and wa as a fraction of the total mass Qy and Qe respectively are identical 


but the absolute mass is allowed to vary. 


subdomain by the total mass and keeping it constant in time [7]. This principle is 
demonstrated in Figure 3.2 for the subdomain w(t) of Q(t). Mathematically the prin- 
ciple may be written as 

1 

al hdx = p(w), (aia) 

6 

w(t) 
where p is a constant for each subdomain, independent of time. Note that in the special 
case where w = Q, then yp = 1. Eq. (3.3) applies a constraint to the PDE (Eq. (3.1)), 
from which the motion of the domain Q and its interior may be inferred. 
To generate the velocity for the interior nodes based upon the relative conservation 

principle, we rearrange Eq. (3.3) and differentiate with respect to time. This can then 


be written in the modified Lagrangian form 


“ ( i » hb) = by(w), (3.4) 


where 6 is the temporal derivative of the total mass 0. The Reynolds Transport Theorem 
[78; 110] (see Appendix B.4) can then be applied to transform it into an Eulerian form 
which explicitly contains the velocity. This velocity represents the deformation of the 


domain and is not necessarily the same as the flow velocity, e.g. u in Eq. (2.26). Under 
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this transformation Eq. (3.4) becomes 


/ hidx + ¢ hv ndD = 6p(w), (3.5) 
w(t) Ow(t) 


where v.n is the normal boundary velocity for the boundary Ow of the region w. By 


applying the Divergence Theorem to the surface integral, Eq. (3.5) can be written as 


i) , {hy + V.(hv)} dx = Op(w), (3.6) 


where v is any sufficiently smooth velocity field that coincides with the normal velocity 


at the boundary. Note that h; can be taken from the original PDE, Eq. (3.1), to leave 
i V.(hv)dx = 6 — i (F(h) +m) dx. (3.7) 
w(t) w(t) 


The rate of change in total mass 6 can be found using the special case w = 2, such that 


Eq. (3.7) simplifies and rearranges to 


oz [ (Fly +m) dt | Vella (3.8) 


= | (F(h) + m) dx + (hv).ndI, (3.9) 
Q(t) AQX(t) 


where we have applied the Divergence Theorem (see Appendix B.5). The rate of change 
in total mass can then be explicitly found provided that h = 0 or v.n is prescribed on 
the boundary. 

In problems of one spatial dimension Eq. (3.7) reduces to an explicit equation for 


the velocity field, 
\ x(t) 
v= fn - | (F(h) +m) to : (3.10) 
0 


for some subdomain (0, 7(t)) of the domain [0, b(t)]. 
In higher spatial dimensions Eq. (3.7) is not sufficient to determine a unique solution 
v, for which an additional constraint is required. Here we follow [21] in prescribing the 
vorticity, the curl of the velocity field v. In particular, by assuming that the flow 
is irrotational, the curl of the velocity field is zero. Therefore there exists a velocity 
potential ~ such that 
v=Vv. (3.11) 


CHAPTER 3. NUMERICAL MODELLING 25 


Substitution into Eq. (3.7) then gives a unique solution for the velocity potential (to 
the extent of a constant), with h = 0 or v.n prescribed on the boundary. The velocity 
field v can then be recovered from Eq. (3.11) [21]. The evolution of mesh nodes x of 


the domain can be found from 
dx s 


ae 12 

ap (3.12) 
Additionally the evolution of the total mass can be determined by integrating 

dog 

—=08 3.13 

dt y) ( ) 


substituting for @ from Eq. (3.9). 
Finally the solution h can be recovered from the relative conservation principle, 


Eq. (3.3), at the new time in the form 


E | has = p(w) = E | ha <J (3.14) 


where t = 0 represents the initial time. The RHS of Eq. (3.14) is known from initial 
conditions. 
We call this the conservation of mass fraction (CMF) moving mesh method and an 


algorithm for this method based upon the above argument is as follows: 


CMF Moving Mesh Procedure 


1. Initially: Compute the total mass 0 (from Eq. (3.2)) and the relative conservation 


constants ju(w) (in Eq. (3.3)) at the initial time. 


2. Calculate the temporal derivative of the total mass 6 from Eq. (3.9), assuming 


that h = 0 or v.n is given on O02. 


3. Find the velocity field v from Eq. (3.7), using a velocity potential in more than 


one dimension. 


4. Advance the nodes of the mesh in time using Eq. (3.12) and a suitable time- 


stepping scheme. 
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5. Update the total mass at the new time from Eq. (3.13), again with a suitable 


time-stepping scheme. 


6. Recover the solution h at the new time step from the relative conservation prin- 


ciple Eq. (3.3). 


7. Repeat steps 2 — 6 until the desired final time. 
This algorithm will be applied to the shallow ice equations in later chapters of this 
thesis. While there is no formal proof available, experience has shown that the mesh 


points tend to move to areas where higher resolution is required, most likely because it 


relies on physical properties of the model to generate the mesh velocity. 


3.4 Adaptive Ice Sheet Models 


In general the use of adaptive methods in the glacial modelling community has been 
limited. The majority of numerical solutions to the ice sheet equations are fixed grid 
methods with either fixed resolution or a predetermined variable resolution. Variable 
resolution methods are useful when the areas that require greater accuracy are known 
in advance; however, the problem of modelling the moving terminus of the glacier is 
not conducive to variable resolution methods. 

The nature of the ice sheet equations can cause the modelling results to suffer from 
numerical problems. The spacing between grid points is critically important, especially 
in certain areas such as the moving front which exhibit steep gradients. However, limi- 
tations in computational power provide a limit on the resolution of the grid. Numerical 
errors may initially start small, but increase non-linearly over time with the high level 
of non-linearity in the mathematical equations [101]. 

There is therefore a strong argument for the use of adaptive grids within dynamical 


ice flow modelling. We summarise the situation below. 


Fixed Grids 


Of the fixed grid approaches there exists multiple ways to discretise the domain to 
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obtain a solution. The most popular is a finite difference type approach, which uses 
a regular grid to allow easy calculation of the ice flux and velocities at the nodes (for 
example {19; 45]). Furthermore, solving the equations on staggered grids, where some 
variables are approximated between the nodes can help to reduce numerical inconsis- 
tencies (see e.g. [86]). 

A small number of methods use a finite element procedure, where the problem is 
solved using a set of equations derived from the weak formulation of the governing equa- 
tions. The finite element method allows for a greater flexibility in the domain structure; 
however, the resulting solutions can be computationally expensive. The best example 
of a fixed finite element method is the model of Antarctica shown in [35; 36; 81], while 
there is also the open source code Elmer/Ice which contains solvers for various problems 
[1]. A finite volume method has also been employed on non-regular grids, such as in 


[53; 60]. 


Adaptive Meshes 

The most common type of adaptivity used are h-refinement techniques. Within 
the ice sheet modelling community there have been a few attempts at applying these 
methods. The first notable attempt appears in 1996 when Lam and Dowdeswell tested 
a number of existing h-refinement methods within a finite volume shallow ice model 
[62]. As in this thesis, the authors were concerned with the accuracy achieved in mod- 
elling the moving front of the glacier. They concluded that while the adaptivity has 
little effect on the steady state solution, the motion of the boundary was significantly 
smoother with less spurious fluctuations than with a comparable unrefined fixed grid 
method. This came at a computational cost roughly 1.8 times higher. 

The first known h-refinement method to include two horizontal dimensions appears 
in Starr [97], where the author demonstrates the potential accuracy gains as well as 
the improvement in computational cost to meet similar levels of accuracy. In addi- 
tion, many h-refinement methods use the AMR (adaptive mesh refinement) technique, 
which contains a fine mesh where required with a series of *blocks’ of decreasing levels 


of refinement moving further away from the dense regions. For example [23] uses the 
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method within a finite element framework to develop error estimates and the likes of 
[39; 24] use AMR to estimate the point at which an ice sheet begins to float, known as 
the grounding line. While not developed specifically for ice sheet modelling, the Parallel 
Hierarchical Adaptive MultiLevel Project (PHAML) has recently been integrated into 
the community ice sheet model (CISM) to allow a combined hp-adaptive method for 
this next generation ice sheet model [112]. 

A location based r-refinement method is introduced in [43], where a one-dimensional 
flowline model is approximated on a stretchable grid. The stretched grid is transformed 
onto a fixed computational domain € € [0, 1], similar to Figure 3.1. The one-dimensional 
version of the shallow ice equations (Eq. (2.14)) can then be written in the transformed 


co-ordinates as 


dh 1 Ohu , b(t) Ah 
b(t) ae | 


where b(t) represents the moving glacial front. A good comparison between this method 


Ot —-(t)-: OE 2) 


and similar fixed grid approaches for approximating the glacial grounding line can be 
found in Vieli [105], which concluded that the stretched grid approach is more robust 
than the fixed grid methods which were sensitive to the horizontal grid spacing. There 
is further development of this stretched grid method applied to the grounding line 
problem in [93]. 

Additionally the Ice Sheet System Model (ISSM) developed by NASA is a location 
based adaptive method that uses the surface velocity to reduce errors [63]. 

At present it appears that there has been no attempt at applying a velocity-based r- 
refinement method to the dynamical ice flow problem. These methods are more widely 
used in the field of oceanography and efforts have been made in the related field of sea 
ice dynamics [107]. 

This thesis will develop the velocity-based CMF moving mesh method introduced 
in Section 3.3 to the shallow ice equations to assess the ability of this type of method 


to simulate dynamical ice flow. 
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3.5 Terminology 


In much of the literature the definition of the dimensions of the model depends on the 
author. For the purpose of clarity we define here the definitions used throughout this 


thesis: 


One-Dimensional Model:- A model that contains one horizontal dimension and is 
vertically integrated to a single layer. In some pieces of literature this is referred 


toasa 13D model. 


Two-Dimensional Model:- A scenario with two horizontal dimensions, again verti- 


cally integrated to a single layer. Sometimes referred to as a 25D model. 


In addition the following definitions are proposed to distinguish between fixed and adap- 


tive methods: 


Grid:- The term ’grid’ refers to the numerical domain for fixed grid methods. The 
structure of the domain is rigid in nature and the grid points are generally evenly 


distributed. 


Mesh:- A ’mesh’ is the numerical domain within adaptive models. It implies the struc- 
ture is more flexible and subject to change, either through additional nodal points 


or the movement of the numerical domain. 
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3.6 Summary 


In this chapter we introduced the concept of adaptive mesh techniques and the scenarios 
in which they provide improvements in terms of accuracy and computational cost, 
compared with traditional fixed grid methods. In particular we focused on r-refinement 
methods, referred to hereafter as moving mesh methods. We defined in detail a velocity 
based moving mesh method that conserves relative mass that will be applied to the 
shallow ice approximation equations in chapter 5. 

We then considered numerical solutions of ice sheet models, examining the previous 
attempts to apply adaptive mesh techniques to the equations governing dynamical ice 


flow. 


Chapter 4 


Data Assimilation 


In this chapter we introduce the concept of data assimilation, a procedure which uses 
a combination of a numerical model and observations to produce a statistically best 
representation of the model variables. We look at the basic aims of data assimilation 
and the initial set up, before discussing three methods that are commonly used in 
practice with applications such as numerical weather prediction and oceanography. We 
then look at the use of data assimilation within the field of glaciology, with a thought 
to the potential observations that are available for use within the assimilation. Lastly 


we look at how data assimilation has previously been used with adaptive mesh models. 


4.1 What is Data Assimilation? 


Numerical models are frequently used in physical systems to provide forecasts of the 
state of the system in the future. While these models can provide good estimates, 
uncertainties in the initial input data often lead to errors as the simulation evolves. 

Similarly it is rare for a complete set of information defining all the required variables 
and parameters to be observed at a specific time. This information, in the form of 
measured observations, also contains uncertainty and random noise in addition to being 
incomplete. 


The idea therefore, is to construct a method to combine numerical models with 


dl 
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observations to produce more accurate representations of the current state of the system. 
This is known as data assimilation (DA) within the fields of climate, atmospheric and 
oceanic modelling. It does however bear close resemblance to the fields of inverse 
problems and control theory. 

Data assimilation methods combine observations with a discrete dynamical model 
of the form 


Zk+1 = Mazz + (Em) k k= 0, 1, ites (4.1) 


where the vector z, € IR? is known as the state vector and contains every variable within 


the dynamical model at the time t,. The model operator M is a nonlinear function 
that defines the evolution of the state variables from time t, to t,4,. The vector (€,,); 
represents the error introduced through the model operator M. If we assume the model 
is perfect [80], then this error may be ignored, with (€,,), = 0. 
In addition, suppose at time t;, there are p observations available that can be related 
to the state vector by 
Vr = Cezy + (€o)k AL o8 (4.2) 


where the vector y;, € R? contains all the observations at time t,. The vector (€,), € R? 


contains the observation errors, which are assumed to be uncorrelated and unbiased [65]. 
It is worth noting that in general the problem is under-determined as there are fewer 


observations available than elements in the state vector, such that p < q. This intro- 


duces the need for the operator C, : R? — R?’, a function known as the observation 


operator which maps the model variables in the state vector to the predicted values of 
the observations. 

For data assimilation the observations y; are statistically combined with our esti- 
mated (or forecast) state vector z/ to gain the most accurate description of the state 


of the system, called the analysis 2°. 


CHAPTER 4. DATA ASSIMILATION 33 


4.2 Data Assimilation Algorithms 


There are numerous variations of data assimilation, each with their own benefits. In 
general these can be classed into two types of method; sequential schemes, which evolve 
the dynamical model (Eq. (4.1)) to a point in time with available observations, whereby 
the assimilation procedure occurs before the model continues, and variational schemes 
which incorporate observations from a window of time simultaneously to improve the 
initial conditions for the dynamical model. 

Here we introduce three methods which are commonly used for operational purposes; 
the first of which is the primary method used in this thesis, while the other two are 
introduced for discussion purposes. The notation used throughout follows [50] where 
possible; however, some symbols are changed to avoid confusion with those used in 


glaciology. 


4.2.1 3D-VAR 


The 3D-VAR method [25; 69; 79] is a variational method based upon a maximum a 
posteriori probability estimate that considers only the spatial dimensions of the system. 
The removal of the temporal variation means the 3D-VAR approach is essentially a 
sequential scheme. The analysis state is found by minimising a cost function that 
balances the model state z; with the forecast state zt and observations y;. This cost 


function takes the form 
1 1 
J(Zx) = 5 (2k — 2} )? BuO! (zy — 2h) + 5 WE — CyZz) Ry" (yn — CeZn) (4.3) 


where the forecast state comes from an evolution of the dynamical model in Eq. (4.1), 


applied to the previous analysis solution. i.e. 


zi = Mze_,. (4.4) 


The matrices By, € R?*? and Rz € R®*? are symmetric, positive definite covariance 


matrices that statistically represent the variances and covariances of the errors in the 


forecast and observations respectively. In essence these matrices control the balance 
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between the forecast and observations and can be thought of as a weighting between 
the two. 
The cost function Eq. (4.3) can also be derived from Bayes’ Theorem; 


p(y |2)p(z) 
plzly) = py © (4.5) 
where the analysis solution is the solution with the maximum probability given the ob- 
servations y and an initial guess z/ [59]. Assuming each of the probabilities in Eq. (4.5) 

is Gaussian, we obtain the cost function in Eq. (4.3). 
To minimise this cost function we use iterative techniques such as the Newton al- 
gorithm, steepest descent or conjugate gradient, which can be computationally costly. 


A control variable transform can be applied to increase the efficiency of the iterations 


and speed up convergence. 


Best Linear Unbiased Estimate 


An alternative method to minimising the cost function is to derive an explicit formula 
for the analysis solution. Differentiating Eq. (4.3) with respect to the model state yields 
the gradient, given by 


VJ (Zx) = By (Ze = zi) + CLR; (ye = CpZr) (4.6) 


where C;, € R’*? is the linearised version of observation operator C;. The analysis 
solution therefore satisfies 


VJ (zg) = 0. (4.7) 


Rearranging Eqs. (4.6) and (4.7) and applying matrix identities leads to the Best Linear 
Unbiased Estimate (BLUE) [65] for the analysis solution 


ay = 2% + Kx(ye — Cxzi) (4.8) 


where K € R?*? is known as the gain matrix [79], given by 


K, = By Cy’ (Cx B,C, + Ry). (4.9) 
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Figure 4.1: Concept cartoon of the optimal 83D-VAR data assimilation scheme, using each 


set of observations sequentially. 


This gives us a solution that is equivalent to iteratively minimising the cost function and 
enables the 3D-VAR problem to be solved sequentially [68]. This is shown pictorially in 
Figure 4.1, where the analysis solution is evolved using the dynamical model in Eq. (4.4) 
to provide the forecast at the next assimilation time. 


We now provide more information about the individual components of Eq. (4.8). 


Observation Operator 


Observations of variables in the system are generally unevenly distributed through 
both space and time and contain gaps where no data is available. It is unlikely that a 
perfect one-to-one correspondence exists between observations y and the state forecast 
z/. The observation operator maps the state vector to observation space so it is directly 


comparable to the observations by: 
Yr = Ch (Zx), (4.10) 


with 


C, : R! > RY. (4.11) 


Observations occur either as a direct measurement of the dynamical model variables or 


indirectly via a mathematical relationship. 
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Using direct observations the observation operator is a linear interpolation matrix 
of size p X q, where each row represents the interpolation of the state vector to an 
individual observation. The elements of each row then correspond to the coefficients of 
the interpolation. 

For indirect observations the operator C;, maps the model variables via a physical 
relationship to be comparable with the observations and is generally non-linear. For 
computational purposes this is generally linearised for use in a data assimilation scheme, 


potentially introducing further error [95]. 


Observation Error Covariance Matrix 


The symmetric, positive definite observation error covariance matrix R statistically 
describes the uncertainty in the errors contained in the observation vector y and is 
defined as 

R, = Elece! |x, (4.12) 


where E|-] denotes the expected value and (€5)% = ¥x —Crzz. These errors arise from a 
number of sources, such as inaccuracies of the measuring instrument, physical features 
that are too small to represent in the model or errors that occur in the observation 


operator. 


Background Error Covariance Matrix 


The matrix B is symmetric, positive definite and referred to as the background error 


covariance matrix. It is statistically defined as 
B, = Eleses |e, (4.13) 


where (€)), = zi —z,, is the difference between the model prediction and the state of the 
system. These errors arise primarily from errors in the initial conditions of the system. 
Arguably the most important part of any data assimilation problem [11], the back- 


ground error covariance matrix By, is also one of the most difficult to determine since 
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these errors are unknown and require approximation. In addition to representing un- 
certainty in the errors in the dynamical model, the elements in B; act to spread the 
information from observations; i.e. how an observation at a point affects the nearby 
points or how observations of one type of variable influence other variables [16]. This 
spreading can be seen in Eq. (4.9) as the matrix B is the last operation in the calcula- 
tion of the gain matrix and is the reason why the background error covariance matrix 


is so critical to the data assimilation scheme. 


4.2.2. Kalman Filter 


An alternative approach is a sequential method known as the Kalman Filter [57; 58]. 
The difference between the Kalman Filter and 3D-VAR is the evolution of the back- 
ground error covariance matrix B. With the Kalman Filter this matrix evolves using 
information from the current analysis and the dynamical model to forecast the covari- 
ances to the next assimilation time. 

The sequential equation for calculating the analysis state vector is the same as the 


one used in 3D-VAR, Eq. (4.8): 


zn = 2h + Ki(yeCnzt), (4.14) 
subject to 
zi = Mazi ,. (4.15) 
The gain matrix is given by 
K, = Pic? (C,P{CT + R,)71. (4.16) 


Here the background error covariance matrix is expressed by P;, to distinguish the 
method from 3D-VAR. After calculating the analysis state vector, the analysis back- 


ground error covariance matrix is calculated using the following equation: 
P? = (I— K,C,)P/ (4.17) 
This is then explicitly forecast forward using the linearised version of the dynamical 


model Eq. (4.4); 
Pi, = MPM’, (4.18) 
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Figure 4.2: Concept cartoon of the Kalman Filter data assimilation scheme, a sequential 
scheme where the background error covariance matrix (dotted lines) is evolved through time 


in addition to the solution. 


and the process begins again. This approach means the error covariances are dependent 
on the model flow and can carry the statistical information about the forecast errors 
through time rather than assuming fixed error covariances (see Figure 4.2). Information 
from all the previous assimilation steps are taken into account, making the Kalman 
Filter an accurate method that is the optimal sequential solution. However in practice 
the CPU power required to forecast the background error covariance matrix means it 


is seldom used operationally. 


4.2.3 4D-VAR 


One of the drawbacks with the 3D-VAR approach is that observations are assimilated 
at set moments in time, which means that unless they are taken at the same moment 
they need to be interpolated to the required time, which results in a temporal error 
being introduced. It is generally computationally unfeasible for sequential schemes 
to assimilate at every time observations occur so we turn to variational schemes to 
incorporate their distribution in time. 

The 4D-VAR data assimilation scheme improves the estimate of the state vector at 


the initial time [64] so that when the dynamical model evolves the state vector better 
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Figure 4.3: Concept cartoon of the 4D-VAR data assimilation scheme, which utilises all 


observations simultaneously to provide the best estimate of the system state. 


matches the available data (see Figure 4.3). Observations can then taken at the time 
they occur until the end of a specified time window, the length of which varies depending 


on the application. The cost function to minimise is then expressed by 


1 _ ly _ 
J(Zo) = 5 (20 — 25)? Bo "(20 — 2)) + 5 Soi — Ci2;) "Rj "(yi — Citi), (4.19) 
i=0 
subject to 
Zi = M(to, ti Zo), (4.20) 


i.e. the state vector at time t; is the solution to the dynamical model evolution of z 
at the initial time t = 0. Eq. (4.19) differs from the 3D-VAR cost function Eq. (4.3) 
by the introduction of a summation over the observations distributed in time. In order 
to minimise Eq. (4.19) with respect to the initial state and thus implement a 4D-VAR 
scheme, the gradient of the cost function requires an adjoint model [31], which steps 
backward in time but involves linearising the dynamical model. In many cases the 
adjoint is difficult to calclulate, particularly when the dynamical model is highly non- 
linear. The solution to the 4D-VAR cost function implicitly evolves the background 
error covariances over the time window. This means that both the Kalman Filter and 
4D-VAR methods give an identical analysis at the end of a time window when the 


dynamical model is linear and the same data is used. 
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4.3 Data Assimilation in Glaciology 


Assimilation in the field of glaciology is still very much in its infancy. While its use is 
common place in atmospheric and oceanic modelling, glaciologists have only recently 
began to look at incorporating observations into their models. 

For example, in [3] observations of the rate of thickness change are used to enforce 
covariance relationships between the rate of thickness change, the flux divergence and 
the accumulation rate. In [106] observations of surface velocity are used to reconstruct 
the flow and stress fields of an Antarctic ice shelf, with the eventual aim of predicting 
the potential collapse of the shelf. 

More recently an adjoint 4D-VAR. style method for a flowline model has been de- 
veloped in [41], where the authors discuss the difficulties in reconstructing the basal 
stresses within glaciers. Another attempt at constructing a 4D-VAR adjoint method 
focuses on trying to find the best estimate of the climatic temperature and volume of 


a flowline model with the aim of extending to a more advanced model in the future [15]. 


4.3.1 Available Observations 


In recent times the observational network of ice has improved dramatically. The use of 
satellites in particular has provided a wealth of information to aid the understanding 
of glaciological processes [4]. 

For example, satellite radar altimetry is able to accurately measure the surface 
elevation of ice sheets [9; 67], while visible imagery can be used to locate the ice extent 
or key features [52]. More recently the CRYOSAT missions are returning accurate 
information regarding changes in ice thickness for land ice and values for the thickness 
of sea ice [30]. The most common observation available is the surface velocity, which 
can be determined by analysing a succession of radar images, although the change in 
velocity with depth will then need to be inferred [54]. 


Sensors fitted to aircraft can provide more detailed information such as the thickness 
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of the ice [70]. In situ ground measurements tend to give very detailed information but 
are generally sparse, both temporally and spatially. In particular the accumulation rate 
of snow to ice is measure by ice-cores and snowpits [103], with some form of interpolation 
required to cover the entire domain. 

Historical data is even more limited, primarily due to ice sheets existing in remote 
areas of the planet. Sparse information of past surface velocity, changes in thickness 


and accumulation rate exist, though it is both spatially and temporally limited. 


4.4 Data Assimilation on Adaptive Meshes 


The use of adaptive meshing techniques within data assimilation methods is also an area 
of interest in its infancy. The application of h-adaptivity and nested mesh numerical 
models can be achieved relatively easily [29], whereas r-adaptive schemes prove more 
difficult. 

It has been determined that location based r-adaptive methods can be incorporated 
into an adjoint solution of a 4D-VAR method with benefits to the overall CPU time 
[33; 34]. The UK Met Office is currently looking at allowing the use of static r-adaptive 
methods within their existing 4D-VAR framework. This is achieved by including an 
additional transformation to simplify the background elements of the cost function 
[88; 89]. 

There is however a lack of study into the use of velocity based r-adaptive meshes 
within data assimilation. These moving mesh methods directly solve the evolution of the 
domain, which we may observe through the use of visible imagery yet this information 


is currently underused. 


4.5 Summary 


In this chapter we introduced the general dynamical model equations and the data 


assimilation problem used within geophysical fields such as meteorology and oceanog- 
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raphy. We have given a brief introduction to three commonly used methods, where we 
discussed the benefits of each. In chapter 6 we will utilize the first of these, the 3D-VAR 
method, to a moving mesh ice flow model. 

We examined the current state of data assimilation in the cryospheric fields and the 
attempts to apply the techniques to dynamical ice flow models. We discussed the avail- 
ability of observations, before detailing the application of data assimilation to adaptive 


mesh methods. 


Chapter 5 


A 1D Moving Glacier Model 


In this chapter we begin by looking at the shallow-ice PDE for modelling dynamical 
glacier flow. We look at how the model is set up and the property of total mass and 
its evolution over time. We then consider the diffusive element of velocity at the front 
of the glacier and the condition upon which this value is non-zero and the scenarios 
to which it applies. We also describe a steady state solution for time-independent 
accumulation/ablation rates. 

Next we describe the application of the BHJ method described in Chapter 3 to 
the shallow ice model from Chapter 2. We detail the extraction of a net velocity 
from the conservation of mass fraction, which results in movement based upon both the 
diffusive flow velocity and the accumulation/ablation terms. From there we will develop 
a numerical approximation to the equations in this moving mesh method using a finite 
difference scheme, before demostrating the method using a simple test scenario. Finally 


we compare the method to existing numerical models using a standardised scenario. 


5.1 The Shallow-Ice Glacier PDE 


For many glaciers, such as Nigardsbreen in Norway [84], the flow is dominated by 
movement along the centre of the glaciers, known as the flowline. It is therefore practical 


to represent these using a one-dimensional model along this flow line. To do this the x- 
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component is set to follow the flow line, starting at the top of the glacier with movement 


in the y-direction ignored. Recall from Eq. (2.14) the mass balance equation; 


m + V.(hu) =m, (5.1) 
which in one-dimension can be simplified to [48]: 

Oh O(hu) 

mie =m. 5.2 

a ae me 


As before, h = h(z,t) is the ice thickness and m = m(z, t) is the ice-equivalent accumu- 
lation rate (including the effect of ablation). In Eq. (5.2), u = ch"*!s” is the diffusive 
flow velocity of an isothermal ice sheet from Eq. (2.26), describing the movement in- 
duced by the ice deforming under its own weight. 

A typical domain for a flow line model is shown in Figure 5.1. The left-hand bound- 
ary represents the ice divide, a section of ice which is connected to either a rock face [61] 
or a further expanse of slower moving ice not included in the model domain [111]. Math- 
ematically we treat this section as static and represent this boundary with a no-flux 


condition at the fixed origin x = 0; 
Oh 
Ox 


Conversely, at the right-hand boundary (referred to as the ice margin, or moving front) 


dx 


Sy 
D> lt 


=0. (5.3) 


=0, u = 
x=0 


xz=0 


x=0 


the ice converges to meet the bed [104], but the location of the boundary moves as the 
front advances or retreats. We will consider glaciers where the gradient at the boundary 
is negative, hz < 0, i.e. the moving front does not overturn on itself. We identify this 
moving boundary x = b(t), with the condition 
h all (5.4) 
x=b 
It is common to non-dimensionalise equations such as Eq. (5.2) that contain physical 
quantities. We detail this scaling in Chapter 7 when the model is presented in two- 
dimensional form with the understanding it is also applied to the 1D model. 
We now consider the change in mass over the entire domain. The total ice thickness, 


denoted here by @(t) and defined as 


| © dx = A(t), (5.5) 
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Height 


0 <b(t) > 
Distance from divide along flowline 


Figure 5.1: Glacier domain along flowline from the fixed ice divide at x = 0 to the moving 


front atx =b 


can be differentiated with respect to time to determine the overall change in mass. 


Using Leibnitz’ integral rule (see section B.1) gives 
d b(t) 


ee 5.6 
=a if hae (5.6) 


4) OR dx 
= / ae + ha = 


Under the boundary conditions (Eqs. (5.3) and (5.4)) the second term is zero. Substi- 
tuting Eq. (5.2) into Eq. (5.7) leaves 

. be) 9 b(t) 

= = 3a [huldx + / mdx (5.8) 
b(t) b(t) 


+ mdz 
x=0 0 


b(t) 
(5.7) 


=—hu 


Using the boundary conditions again forces the first term to be zero, since at the ice 
divide we require the ice to remain stationary, implying u = 0. This leaves the change 


in ice thickness over the domain as 


. b(t) 
= mdz. (5.9) 
0 


As a result any change in the total ice thickness over the whole glacier is due solely 
to the source term in Eq. (5.9), ie. the global change in ice thickness equates to the 


net accumulation/ablation over the whole glacier. 
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5.2 Characteristics of the Shallow Ice PDE Model 


We now examine features of the model that provide insight into the behaviour of glaciers 
under certain conditions. Firstly we consider the diffusive flow component and its 
impact on the moving glacier front before considering a steady state solution where the 


ice thickness profile remains constant over time. 


5.2.1 Diffusive Velocity at the Glacier Front 


Glaciers are unique entities as they facilitate an overall advance or retreat of their 
expanse despite the fact that the diffusive flow is unidirectional. Naturally there will 
also occur periods where the front is stationary, either at the cross period between 
advance and retreat or simply a momentary pause in movement (momentary on a 
glacial time scale). Here we assess the movement of the glacier front generated by the 
diffusive flow velocity u in Eq. (5.1). 

Under the Shallow Ice Approximation detailed in Section 2.4.2 the 1D diffusive 


velocity u in Eq. (5.2) can be written as 
Caches. (5.10) 


where c = —2Ap"g"/(n + 2), under the assumption that temperature and density are 
constant throughout the glacier. This means that the parameter c can be written as 
a single negative constant. The surface elevation, s, is the total height of the glacier 
including the topographical bed as well as the ice thickness. For simplicity we begin by 
making the further assumption that the bed is flat, with s = h, so that 


ieee iiaalas ie (541) 


Expressing the diffusive velocity in the form of Eq. (5.11) suggests a problem when 
applying the boundary condition h = 0 at x = b(t), apparently yielding a zero diffusive 
velocity and resulting in a boundary that only moves as a result of accumulation or 


ablation. However, it is perfectly possible for the diffusive velocity u to be non-zero 
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when h = 0 and h, is infinite such that h"t'h” is finite. 
To demonstrate this we shall perform an analysis on the local profile of the flow 
velocity, close to the moving front. We begin by utilising the chain rule of differentiation 


to write Eq. (5.11) in the more useful form 


n n n n ‘ n n ih 
u=c(h™/p,)” =c (= :] ase a) (5.12) 


from which it can be seen that the diffusive velocity is finite at the boundary provided 


that (h@rt+)/”),. is finite. 


We now wish to analyse the conditions under which Eq. (5.12) provides a non-zero 
velocity at the boundary. To achieve this we note that the profile of ice thickness must 
satisfy the boundary condition h = 0 at x = b, and therefore at any moment in time 
it must contain the factor (b — x). We may then let a = a(t) > 0 be the largest real 


number such that h close to the moving boundary b can be written in the form 
h= (b—2)*g(x) (5.13) 


where the function g(x) > 0 is finite and has a finite derivative at the boundary x = b. 
Then 

penth)/n = ((b = x)%g(x))Cnt vin a (b = ag) 222+ / G(x) ; (5.14) 
where G(x) = g(x)@"+)/" > 0 is also finite and has a finite derivative at x = b. The 
value of u in Eq. (5.12) at the boundary is proportional to the derivative of Eq. (5.14), 


which assuming continuity of wu and G in the vicinity of the boundary, can be found by 


taking the limit as x tends to b 


lim(hOnth/"), = lim [{(b — c)°@rt DG (a)} J. (5.15) 


xb x2—b 


We then note that since both 


lim [(b — 2)°P"tY/G(x)] =0 and lim (b— 2x) =0, 


a—b a—b 


we may employ L’Hopital’s rule (see Section B.2), leading Eq. (5.15) to become 


_ p)a(2n+1)/n 
lim(An@rt)/™), so, te lim (b 7 - D G(x) ee lim(b —_ g Cnn TG (a); (5.16) 
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It follows that the limit of the diffusive velocity (Eq. (5.12)) is proportional to 


lim (psy: (5.17) 


ab 


Since n is a given exponent, there exists a critical value for the power a, denoted 


Q@-, which leads to the following properties for the diffusive velocity u at the boundary 


as x > b; 
Casel: a>ae => u->Oaszr—>bd, (5.18) 
Case 2: a@=a,. =>  wremains finite as 7 — b, (5.19) 
Case 3: a<a, => wu—>coasxz—b (momentarily). (5.20) 


where the critical value of a is given by 


c= Bae): (5.21) 


Thus the ice diffusion induces no movement at the boundary when the ice thickness 
profile given in Eq. (5.13) has a value of a in the range defined by Case 1. As the 
profile changes over time and a tends to the critical value a, in Eq. (5.21), the Case 
2 condition is met and the diffusive velocity is then non-zero. If the thickness profile 


changes to the form of Case 3 then the velocity is momentarily unphysical. 


Specific Value for the ice flow exponent, n 

At this point we focus on a specific value for the ice flow exponent, n = 3, which is 
an empirically found standard value in ice sheet modelling [102]. Eqs. (5.18) to (5.20) 
can then be rewritten using Eq. (5.21): 


3 

Case 1: eas u>Oasz—>b, (5.22) 
3 

Case2: a= go Cee al remains finite as x > b, (5.23) 
3 

Case 3: a< 7 => wu-—>ooasx—b (momentarily). (5.24) 


We may demonstrate the impact of Eq. (5.21) with n = 3 by using an initial thickness 
profile with the simple function g(x) = (1+ 2)® in Eq. (5.13) with the boundary b = 1, 
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such that 
h=(14+2)*%(1-—2)*. (5.25) 


By varying the value of a we see the diffusive velocity at the boundary up, is discon- 
tinuous at the moment a reaches 3/7 from above (Figure 5.2(a)), jumping from zero 
to a finite value. In Figure 5.2(b) we see u in the vicinity of the boundary for varying 
values of a. As we approach the critical value from above the discontinuity forms until 
a reaches 3/7 itself, at which point a non-zero finite velocity exists and the boundary 


begins to move as a result of diffusion. 


a=3/7 

~ - -a= 49/112 

a = 25/56 

— a= 13/28 
® @et@ [ome ee ee 

+ = 15/28 OE Bi 


b 
Diffusive Velocity 


a , , Distance from Divide 
(a) (b) 


Figure 5.2: Varying the parameter a in Eq. (5.13) with g(x) = (1+ 2x)® leads to a critical 
value when considering the diffusive velocity at the glacier front. a) As a > 3/7 from above 
the boundary diffusive velocity uy jumps to a finite value. b) Velocity in the vicinity of the 
boundary for varying values of a. As a — 3/7 the velocity profile steepens, maintaining a 


zero velocity on the boundary until a reaches 3/7, at which point it jumps to a finite value 


Continuing with this specific example, we consider the diffusive velocity over the 
whole domain (Figure 5.3(a)) where we see that not only does the value of a affect the 
velocity of the boundary, but the velocity profile as a whole. We see the peak value for 
the diffusive velocity moves from close to the centre of the domain when a = 1, to the 
boundary as @ approaches 3/7. When it reaches this value the peak velocity is located 


at the boundary itself. In Figure 5.3(b) we see that when a reaches the critical value 
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Diffusive Velocity 


Ice Thickness 


4 05 06 
Distance from Divide 


(b) 


Figure 5.3: (a) Diffusive velocity. As a — 3/7, the peak velocity approaches the boundary, 
and when a = 3/7 (red line) the boundary velocity is finite whereas all other values have a 
zero boundary velocity. (b) Ice thickness profiles for a = 1 (solid) and a = 3/7 (dashed). 


Notice the infinite gradient at the boundary when a = 3/7 


the gradient of the ice at the boundary b = 1 is unbounded, which we highlighted as a 


requirement to get a non-zero velocity in Eq. (5.11). 


Incorporating the glacier bed 
The above analysis has been performed under the assumption that the bed of the 
glacier is flat. In reality this is rarely the case, so we may define the surface elevation 


as a summation of the bed and ice components 
Sa hy (5.26) 


where z = z(x) represents the topography under the glacier. Continuing with the 
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assumption that n = 3 allows us to write the diffusive velocity from Eq. (5.10) as 


u = ch*[(z + A)z]* (5.27) 
= ch’ |Z, + h,|° (5.28) 
= cz3h* + 3c2z2h*h, + 3cz,h*h? + ch*h?. (5.29) 


Note that the last term in Eq. (5.29) is the same as the single term in the flat bed 
scenario, Eq. (5.11). Expressing the terms in Eq. (5.29) in the same form used in 
Eq. (5.12) we get 


9c 


3c c 
3,4 2/5 = 3) \2 
u=czh*+ 7 z(h?)e + grei(h er 313 


COR ake (5.30) 


The inclusion of a bed does not change the structure of the ice thickness, so again 


we may substitute in the expression for h given in Eq. (5.13) 


w= cah(b— 2)" 9(o)' + F2X((b- 2)" 4(2))e 
+ Szaf (b= 2)*0(e))eP + gee ((6— 2)" 9(a)™),}. (6.31) 


The diffusive velocity at the moving boundary is found by taking the limit of wu as 


x— b, 
tim u = lim |c28(b— 2)g(2)* + £22((b — 2) *9(2)")e 
+ Szef (b= 2)*9(2))eF + aot (( = 2)g(2)")_}5| - (6.32) 


The first term tends to zero at the boundary and again we apply L’Hopital’s rule to 

the remaining three terms individually to leave 

BE! aes mea! c : —_ 
u(b) = —= 2 lim|(b — x)" g(x)"] + s20[— lim{(b — 2)" g(x) "}]° 
b 3 rb 
9 
+ sra[- lim{ ((b- 2)"°-1g(2)")}. 6.33) 
— 

The analysis of the non-flat bed scenario leads to the same critical value of a defined 


in Eq. (5.21), as the following argument shows. Using the same technique here, but 


treating these terms separately leads to a situation where each term provides a different 
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critical value of a to the one found previously. The first term yields the value of 1/5, 
the second 1/3 and since the last term is the same as the flat bed, its critical value is 
3/7 as before. 

Since a cannot go below any of these values without encountering an infinite velocity 
(see Eq. (5.20)), the lower limit for a must be the largest value, namely 3/7. There- 
fore the flat bed analysis holds and the same condition for a diffusive velocity on the 
boundary exists regardless of whether there is a topographic bed present. 

In addition, while we made the assumption in Eq. (5.10) that the temperature and 
density were constant, the analysis above is in fact independent of these parameters 
since the analsis applies to a factor of h. In other words, if temperature and density 
varied spatially they could be included within g(x) which does not impact the condi- 
tion upon a. Differing these values will provide a different value for the finite velocity 
achieved when Case 2 is met; however, the boundary velocity will remain zero while in 
Case 1 regardless of the temperature or density. 

We now have a condition upon the ice thickness profile which allows the diffusive 
flow of the boundary to be non-zero. In the absence of any accumulation/ablation this 
is the sole cause for movement of the ice and it would be possible to determine when 


the glacier remains stationary and when movement will begin. 


5.2.2 Steady-state solution 


It is possible for glaciers to enter an equilibrium state where the net mass of ice remains 
constant over a period of time. We now introduce a solution to this equilibrium state, 
known as a steady state solution, which will be used to test the accuracy of our numerical 
model later in the chapter. 

One approach to finding this steady state solution is to rewrite the problem as an 
energy functional and minimise, see [56], though here we will directly find a solution 
by integration. Glaciers enter steady state when the amount of ice accumulated is 


balanced by the amount lost through ablation. Mathematically we require the change 
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in ice thickness to be zero, h; = 0, which reduces Eq. (5.2) to 
(h(x )u(x))2 = m(z), (5.34) 


where we assume the ice equivalent accumulation rate m to be independent of time in 


order to maintain equilibrium. By integration we may write Eq. (5.34) as 


ies [ m(y)dy (5.35) 


where y is a dummy variable. Substituting the flat bed diffusive velocity wu from 
Eq. (5.11) gives 
cans = f m(y)dy, (5.36) 
0 


which we may write in an alternate form similar to Eq. (5.12) 


c ( we 3) {(n@rt2)/n), y= fo my (5.37) 


Next we move all the constants over to the right hand side before raising both sides to 


(Henin), = (=) ( i, “n(y)dy) ' (5.38) 


Integrating a second time gives 


pentay/n (2) [ ( [ome mi yee)” dy (5.39) 


where z is a second dummy variable of integration. Finally raising both sides to the 


the power of 1/n, 


power of n/(2n + 2) gives the ice thickness at any point x in steady state 


ney (BVT F (fmerae) a]. a 


In order to reconstruct the ice thickness profile using Eq. (5.40), we require the 
location of the moving boundary b, which in steady state will be stationary. To find 
this value we require the net accumulation and ablation across the whole domain to 
be zero, so that although the glacier ice is still flowing the overall movement is zero. 


Therefore in steady state: 


b 
) mdz = 0 (5.41) 
0 
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With Eqs. (5.40) and (5.41), if we know the ice equivalent accumulation rate across 
the entire domain or we choose a function to represent it that is time independent, 
we can find the location of the boundary and thus construct the corresponding steady 
state profile. We make use of this profile to test the results of our moving mesh scheme 


in Section 5.5.1. 


5.3 Glacier Movement 


We now turn our attention to constructing a method that describes the overall move- 
ment of the glacier. This net velocity incorporates the ice diffusive velocity and accu- 
mulation/ablation effects. The method for finding this velocity is based upon the BHJ 


method discussed in Section 3.3. 


5.3.1 Conserving Mass Fractions (CMF) 


For a time-dependent PDE with a source term, the relative mass of a local subdomain 
is conserved by enforcing the conservation principle (Eq. (3.3)) over time. For glaciers 


with ice thickness h this may be written mathematically as 


a(t) = 
hd i ee 

Io. hdx _ if hdz = (2), (5.42) 
ie hdx 4t) Jo 


applied to any moving subdomain (0, 7(t)) of [0, b(¢)], with 0 < & < b. The variable 
6 in Eq. (5.42) represents the total mass defined in Eq. (5.5). We remark that when 
described in this manner the constant-in-time function u(Z) € [0, 1] forms a cumulative 
function with j.(0) = 0 and yu(b) = 1. This method is referred to as the Conservation 
of Mass Fractions (CMF) approach throughout this thesis. 


To extract a velocity we begin by differentiating Eq. (5.42) with respect to time to 
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leave 
a) dO 
— = Ll 5.43 
dt i: aed arr a 
b(t) 
=U if mdz, (5.44) 
0 


where we refer back to the differential property of total mass Eq. (5.9). Carrying out 
the time differentiation on the left hand side of Eq. (5.44) using Leibnitz’ Integral Rule 


qd #0 a(t) it 
_ hide = hid. h— 
didi. | es ( _ 


Applying the boundary condition in Eq. (5.3) forces the velocity dz/dt to be zero at 


gives 
a(t) 


(5.45) 


0 


the ice divide x = 0. Substituting for h, from the shallow ice PDE Eq. (5.2) leaves 


qd peo a(t) ds 
ae ae —(hu)» we 
ds dx i (—(hu)e +m) dx + (: z) 


a(t) 
a(t) ae 
= mdz + (5 = in) (5.46) 
a(t) 
where again the boundary condition in Eq. (5.3) requires the flow velocity u to be zero 


at the fixed boundary « = 0. Combining Eq. (5.44) and Eq. (5.46) we obtain the net 


‘ dt 


velocity of an arbitrary interior point 7(t) to be 


de le Io? mdx — fr mde| 


apo h 


(5.47) 


This velocity describes the net overall movement of a point in the glacial domain as a 
combination of the flow velocity wu and the accumulation rate m. Using this we may 


determine the changes in the domain over time. 


5.3.2 Net Velocity at the Glacier Front 


The velocity given by Eq. (5.47) describes the movement of any interior point. Assuming 
continuity of the velocity the corresponding net velocity at the moving front at x = b 


may be indentified by taking the limit of Eq. (5.47) as Z(t) > b; 


dt tad b(t) a(t) , 
lim — = im — — ‘ A 
a(t)>b dt ul, + a(t)>b h vf nee / race Cee) 
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The limit is undefined since h > 0 as X —> 6b; however, since the quantity in square 
n~ 5 : 5 : : : dx 
brackets also + 0 as x(t) — b, then by l’Hopital’s Rule (assuming continuity of $ and 


m in the vicinity of the boundary); 


bic it : BO E fe max — fF mde or 
a(t) 3b dt uly + a(t) 3b hy oe?) 
Ppa 
= li _L oaeameaescies 
ul, + ek ie (5.50) 
m 
= (45> = (5.51) 


This is the velocity at the boundary, which incorporates the accumulation m (nor- 
malised by the local gradient h,.) as well as the diffusive velocity u. Since h, < 0 at the 
front x = b, the second term in Eq. (5.51) is positive or negative depending on whether 


m is positive or negative. 


If the gradient h, at the front is shallow, the net boundary velocity is mainly in- 
fluenced by the accumulation/ablation component and allows the front to advance or 
retreat. Conversely, for steep gradients the accumulation term becomes less important 
and the flow velocity u dominates. This means that the front is susceptible to the 
analysis of the diffusive velocity in Section 5.2.1: in this case the overall movement of 


the glacier front can then only advance or stay, not retreat. 


5.3.3 Radial Adjustment 


Some glaciological features, such as an Ice Dome (a mass of ice with a rounded, gently 
sloping dome [2]) can crudely be represented using a one dimensional radial model. It 
is also beneficial to create idealistic radially symmetric configurations to obtain steady 
state solutions for higher dimensional models |92]. For a radially symmetric icesheet 


the shallow-ice PDE (Eq. (5.2)) can be adjusted to give 


Oh | 1O(rhu) _ 
dtr or ee?) 
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where the diffusive velocity (Eq. (5.11) becomes 


2Ah nin {Oh\” 


In addition the total mass defined in Eq. (5.5) is now calculated using 


b(t) 
O= i hrdr, (5.54) 
0 


whose rate of change is 


dg Xt (0) 
O= al hrdr = | mrdr. (5.55) 
0 0 


The CMF approach of Eq. (5.42) is based upon the assumption that relative mass is 


conserved. In a radially symmetric environment this implies that 


F(t) 
=f rhdr = v(r), (5.56) 
0 


where the constants v are independent of time and represent the radial equivalent of 
the constants y (Eq. (5.42)). These are defined on any moving subdomain (0,7(t)) of 
(0, b(t)|] leading, in a similar way as in Section 5.3.1, to a similar expression for net 


velocity 
as E i mrdr — io merdr] 
co h 


With this small adjustment we are able to simulate a radially symmetric domain using 


(5.57) 


a one-dimensional approach. 


5.4 Numerical Approximation 


With the method for extracting the net movement of a glacier outlined in Section 5.3 
we now define the techniques to numerically approximate this approach in the form 
of a moving mesh method. We begin by detailing the numerical flowline model before 
examining the numerical approximation to the velocity at the glacier front and finally 


the numerical radial model. 
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5.4.1 A Finite Difference Scheme 


We now describe a finite difference numerical scheme based on the CMF approach in 
Section 5.3. Denoting the mesh nodes by {X;} for (¢ = 1,2,...,.N), the initial domain 


may be discretised such that 
0 = X1(t) < Xo(t) <--- < Xy(t) = W(t). (5.58) 


We define H;, U; and j1; as the discrete approximations of h(X;), u(X;) and w(X;) re- 
spectively. To initialise the moving mesh scheme the domain {X;} and the ice thickness 
{ H;} is required. The numerical scheme forms an algorithm that can be split into three 


stages; 


Stage 1 
The first stage of the algorithm approximates the mesh velocity from Eq. (5.47) using 


finite differences to give a discrete velocity at each interior node in the form 


dX. {1 ~ mdx — Pane mdz} 
v U; 


aM = u,- 7 (5.59) 
nm \" (Hern pentdin\” { ps fo? mda — fe mdz} 
ee he (54) X; — Xj_4 7 Hi; : 
(5.60) 


where the derivative of h@’+)/”" in u is approximated using an upwind difference and 


the integrals are approximated by the composite trapezium rule (see Section B.3). 


Stage 2 
In the second stage we advance each of the nodes X; forward in time using an explicit 


Euler scheme, 


k 
XP aXe oe At (=) ; (5.61) 


for all i. Here k denotes the time discretisation level. The same approach is applied to 


the total mass 0**1. Using Eq. (5.9), the explicit Euler equation for updating the total 
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mass is given by 
pk+1 


okt! — ok + At i? mdz . (5.62) 
0 


where again we approximate the integral using the trapezium rule. 


Explicit Euler time-stepping is sufficient for our purpose provided that the time step 
is sufficiently small to ensure stability and prevent the nodes from crossing each other’s 
path. While there is no formal definition for a moving mesh time-step, we may use as 
a guideline the general rule for diffusion equations: 


(Az)? 
a 


pee (5.63) 


where D = ch”*?|s,|"—" is the diffusivity of the PDE. 


Stage 3 

After finding the new locations for each of the nodes, the final stage of the algorithm 
recovers values for the ice thickness H; at each of the nodes using a discretised form of 
the conservation principle (Eq. (5.42)). Since the integral in Eq. (5.42) is constant in 


time for all z(t), we may approximate this equation in the form 
1 ee k+1 
grti if hae] = fig. — Mi-1 (5.64) 
where the limits in the integral arise from an incremental form of Eq. (5.42). Approxi- 
mating Eq. (5.64) by a midpoint rule, we obtain 


k+1 __ label e Cree = pia) 


See ey 


(5.65) 


which provides the ice thickness at the next time step. Eq. (5.65) forces the relative 


mass to be conserved in time, thus maintaining the initial principle of the method. 


5.4.2 Numerical Diffusive Velocity at the Glacier Front 


We now relate the numerical approximation to the diffusive velocity at the front to 


the analysis in Section 5.2.1. Recall that the upwind approximation U; to the diffusive 
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velocity u(X;) in Eq. (5.60) is 


n (2n+1)/n (2n+1)/n\ ™ 
n H. Si 
U;=- : : 5.66 
Hence the approximation to the diffusive velocity at the glacier front Uy is given by 
4 " (Hy_1)24! 
Ure ‘ 5.67 
a (sn) Crem. eee 


since Hy = 0. At first glance this bears little resemblance to the analytic diffusive 


velocity at the boundary Eq. (5.17), where 
ube baa) ere: (5.68) 


But by rewriting this expression, using Eq. (5.13) we obtain 


(b _ gyrate 


u(b) x (6-5) (5.69) 
_ han +1) 
ae (5.70) 


since h x (b— x)*. Eq. (5.70) has a very similar form to Eq. (5.67). Asymptotically 
when Xy_, is close to Xy (= 6) and a > n/(2n 4+ 1), the expressions in Eqs. (5.67) 
and (5.70) are small, but this is not the case if a = n/(2n + 1) as demonstrated in 
Section 5.2.1. 

Therefore the closer together the nodes are at the glacier front, the more effective 
Eq. (5.67) is at providing a good approximation to Eq. (5.70). We may then say that 


the numerical scheme closely reflects the analytic behaviour of the diffusive velocity. 


5.4.3. Radial Adjustment - Numerics 


The numerical approximation to the radial adjustment made in Section 5.3.3 is very 
similar to the non-radial scenario in Section 5.4.1. The radial domain is discretised into 


a mesh {R;}, where (i = 1,2,...,N) and 


0 = Ri(t) < Rot) <--- < Ry(t) = b(t). (5.71) 
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We define H; and v; to represent h(R;) and v(R;) respectively. Again the algorithm 


can be split into three stages; 


Stage 1 


Stage 1 of the algorithm approximates the velocity Eq. (5.57) by finite differences to 


give 
n b(t) Ri(t) 
dR; - n ee _ BO {v af mrdr — is mrdr } 
dt “\On41 Roh H, 


(5.72) 


at each node R;. 


Stage 2 
The second stage advances the mesh nodes R;(t) with this velocity by the explicit Euler 


scheme . 
dR; 
Re RY At ( ) (5.73) 
dt 
for all i. The total mass ¢ is also updated using an explicit Euler equation on Eq. (5.55) 
to give 
peti 
ort = oF EAT | mrdr . (5.74) 
0 


where we approximate the integral using the composite trapezium rule. 
Stage 3 


Finally the discrete ice thickness H; is recovered algebraically from the incremental 


form of Eq. (5.56) 
1 
perl 


using a midpoint discretisation: 


Rixi(t) pas 
': rhdr = Vix, —Vi-1, (5.75) 
Ri_1(t) 


k+1/,,. BF ae 
HF = ae nai (5.76) 
i+1 aed 4—1 
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As with the non-radial approximation the relative mass is conserved in time and the 


algorithm is repeated at each time step. 


5.5 Experiments 


We now examine two scenarios to evaluate the moving mesh method. The first is a 
simple test scenario designed to demonstrate the method covering the core aspects of 
glacier flow, while the second is a radially symmetric scenario designed to test and 


compare the efficiency of numerical ice sheet models. 


5.5.1 Testing the Model 


For an initial test we shall examine the ice thickness profile encountered during the 


diffusive velocity analysis (Section 5.2.1). The initial profile is of the form 
h = (b* — x) (OshL) 


on an initial domain « € [0,b] where x = 0 is fixed and the initial moving front x, = 1. 


For the source term m we use the quadratic function 
m=y(1- =), (5.78) 


where ( defines where the source term changes from accumulation to ablation (called 
the equilibrium line) and y is a parameter to control the scale of the source term. 
This function is independent of time and as such there exists a steady state solution 
as detailed in Section 5.2.2. Using Eq. (5.41) it can be shown that in steady state the 


boundary, b,,, will be located at 
Dug =AP80: (5.79) 
Using the parameters defined in Table 5.1 we may test the numerical model and its 


behaviour. 
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Physical Parameters 


3 Flow-law exponent 

AS 10 (Pa) a Flow-law parameter 

g=9.81ms~? Acceleration of gravity 

p = 910kgm-? Ice density 

c= —2Ag"p"/(n + 2) Constant Parameter 

y = 0.0005 Scale of accumulation rate 
Computational Data 

IN’ = 26 Number of gridpoints 

Ax’ = 0,02 Initial grid spacing 

AGS 05 Time Step 

T Final Time 


Table 5.1: Values of the physical parameters in the test model and the computational data 


used in the test equations. 


Glacier in Advance: a = 3/7, 6 = s, 2 = 120000 

By choosing the equilibrium line 6 = 1/2 such that the boundary of the steady state 
solution is greater than the initial boundary we can demonstrate an advancing glacier. 
We set a = 3/7 to satisfy the condition on the diffusive velocity for initial movement 
in Eq. (5.21). We observe in Figure 5.4(a) an overall increase in the amount of ice 
in the domain, along with movement of the front towards the steady state boundary 
b., = V1.5 defined in Eq. (5.79). Since there is an initial diffusive velocity the movement 
of the numerical nodes of the mesh are all in the same direction (see Figure 5.5(a)). 

Conversely if we set a = 1 so that the diffusive velocity on the boundary is initially 
zero, we observe that the boundary actually begins to retreat due to ablation from the 
source term Eq. (5.78). However as the ice thickness profile changes over time to satisfy 
the conditon on the diffusive velocity at the boundary, the glacier begins to advance to 


reach its steady state profile (see Figures 5.4(b) and 5.5(b)). 
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Ice Thickness 


02 05 
Distance from Ice Divide 


(a) Advancing Glacier; a = 3/7, 8 = 3 


Ice Thickness 


02 06 
Distance from Ice Divide 


(b) Advancing Glacier; a = 1, 8 = $ 


Ice Thickness 


02 05 
Distance from Ice Divide 


(c) Retreating Glacier; a = 3/7, B = 75 


Ice Thickness 


02 06 
Distance from Ice Divide 


(d) Initially Stationary Glacier; a = 1,6 =1 


Figure 5.4: 1D test equations: Ice thickness over the domain evolving from the initial profile 


(green) to a final steady state profile (blue). 


Glacier in Retreat: a = 3/7, 8 = 4, T = 6000a 


Let us now choose the equilibrium line to be located close to the ice divide (x = 0), 


for example, 3 = 3/10. Under this condition the location of the steady state boundary 
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(Eq. (5.79)) is to the left of the initial boundary. We therefore expect the front of the 
glacier to retreat and the ice thickness to decrease. 

In Figure 5.4(c) we see a rapid initial retreat, where the boundary moves further 
back than the steady state location. The glacier then appears to stabilise and advances 
back up to the steady state location. This is shown in the nodal characteristics plot 
(Figure 5.5(c)) where we observe the initial retreat followed by a slow advance back to 
the steady state boundary. Crucially we see the interior points also retreat in line with 
the boundary which means the nodes avoid crossing each others paths, an important 


stability property for any moving mesh method. 


Stationary Front: a=1,6=1, T = 1500a 


Selecting 8 = 1 in Eq. (5.78) yields zero accumulation/ablation at the initial bound- 
ary b = 1. Under this condition the diffusive velocity becomes the sole contribution 
to the velocity at the boundary, Eq. (5.51). Therefore the boundary responds to the 
analysis in Section 5.2.1 and by setting a = 1 we ensure the boundary is initially sta- 
tionary. 

It is clearly visible in Figure 5.4(d) that the profile evolves towards an infinite gra- 
dient at the front to satisfy the finite diffusive velocity criteria in Eq. (5.19). Once this 
criteria is met the front begins to advance. Examining the movement of the mesh nodes 
in Figure 5.5(d) we can see that the interior nodes move towards the boundary as the 
interior ice is flowing, until this boundary profile is satisfied and the boundary itself 


begins to move. 


Summary 

We have shown that the numerical approximation described in Section 5.4.1 can 
competently model a glacier that is advancing, retreating or remaining stationary. The 
method itself is computationally cheap, each scenario running in under a minute for 
the experimental data given. We now turn our attention to a practical example for 


comparison with existing glacier models. 
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(c) Retreating Glacier; a = 3/7,8 = 4,T = (d) Initially Stationary Glacier; a = 1, 8 =1,T = 
6000a 1500a 


Figure 5.5: Evolution of the mesh nodes over time, including the boundary point (red) 


representing the glacier front. 
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5.5.2 European Ice Sheet Modelling INiTiative 


The European Ice Sheet Modelling INiTiative (EISMINT) experiments [49] give a set 
of standard benchmark scenarios to test and compare numerical models of ice sheet 
behaviour. The initial paper [49] contains information about two scenarios; the first 
is a fixed margin experiment on a square domain which deals with calculating the ice 
discharge at the boundary and the second is a radially symmetric situation with a 
moving boundary. We shall replicate the second experiment using the radially adjusted 
CMF method detailed in Section 5.3.3. 

The experiment begins with a flat bed with zero ice present. The CMF method is 
dependent on preserving relative mass, for which we need an initial mass to calculate. 
To achieve this we take the ice thickness profile after one time-step, which corresponds 


to the source term only as there is no ice to diffuse. This term is given by 
m = min{0.5, y(8 — x)}, (5.80) 
and so the initial profile of ice thickness is 
h° = At x min{0.5, 7(6 — x)} (5.81) 


with x € [0,4.5 x 10°]. In Eq. (5.81) the equilibrium line is located at 6 = 4.5 x 10°m 


1 


and the scale of the source term is given by y = 10-°a~!. Again the source term is 


independent of time, so there exists a steady state solution. Using a radially adjusted 
version of Eq. (5.41) 
b(t) 
| nrdr = (5.82) 
0 


we may calculate the steady state boundary location to be 
bs = 310.81 ki, (5.83) 


and the steady state ice thickness may be reconstructed using the method in Sec- 


tion 5.2.2; 


CHAPTER 5. A 1D MOVING GLACIER MODEL 68 


Physical Parameters 


n=3 

A= 10 °(Pa) a * 
g=9.81ms~? 

p = 910kgm-? 

m = min{0.5, (6 — x)} 


y=10 ma tkn 


Flow-law exponent 

Flow-law parameter 

Acceleration of gravity 

Ice density 

Ice equivalent accumulation/ablation (source) 


Slope of source function 


B = 450km Equilibrium line location 
Computational Data 

N= 16 Number of gridpoints 

Ar® = 3 x 10*m Initial grid spacing 

At=2a Time Step 

S23 <0 Final Time 


Table 5.2: Values of the physical parameters in the model and the computational data used 


in the EISMINT scenario. 


The physical parameters for the model are provided in Table 5.2 along with the numer- 
ical data. For direct comparison with the fixed grid methods grid methods presented 
in [49], performed on a 31 x 31 regular rectangular grid, the 1D moving mesh model 
uses 16 grid points, initially spaced evenly. This allows for a time step of 2a, up to 
a final time of 30,000a. Note that the unit a refers to one annum, equivalent to one 
year. Computationally this solution takes less than five seconds to run on a standard 
desktop. 

As the problem is radially symmetric, the results of the one-dimensional flowline 
method may be presented as a circle which allows for comparison with two-dimensional 
methods. 

The numerical solutions to the experiment in [49] are found using traditional fixed 
grid methods on evenly spaced grids. This means that they require some form of ex- 
trapolation or interpolation to find the boundary location as the boundary generally 


falls between two grid points in a 1D flowline model or four with a 2D grid. As a 
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result fixed grid methods on a regular two-dimensional rectangular grid, such as those 
presented, do not return a perfect circle due to the location of the grid points, as shown 
in Figure 5.6(a). This is not an issue with the flowline models such as the one used 
here, which when rotated naturally give a circle (see Figure 5.6(b)). 

Instead a direct comparison between flowline models would provide more of an in- 
sight. When comparing, we see that the moving mesh method is able to get significantly 
closer to the exact ice thickness profile than the equivalent fixed grid method in Fig- 
ure 5.6(c), especially close to the moving front. We also see from Figure 5.6(d) that 
the diffusive velocity in steady state is similar between the two approaches, with the 
main difference arising at the boundary, where the diffusive velocity can be explicitly 
calculated in the moving mesh approach. The fixed grid schemes require interpolation 
to calculate this value. 

Expressing these results in Table 5.3 shows that the CMF moving mesh solution 
is able to get much closer to the exact boundary position than the average fixed grid 
solution, while the thickness at the ice divide (where r = 0) is slightly higher than 
both the fixed grid and exact solutions. This is most likely due to an underestimation 
of the flow velocity in the accumulation regions of the domain, which is visible when 
compared to the fixed grid in Figure 5.6(d). A slower velocity transports less ice to 
the areas where mass is lost due to ablation resulting in an overestimation of the ice 
thickness in the accumulation regions, where the divide is located. 

We have demonstrated that the moving mesh method provides a better approxima- 
tion than traditional fixed grid methods to a flowline model using the same number of 
grid points. We have also tested the model using scenarios containing time-dependent 
source terms, such as the sinusoidally varying expression found in [49]. In these cases 
we found that the moving mesh method produces a good approximation and is able to 


deal effectively with a domain that is continuously advancing and retreating. 


Convergence 
An important test of any numerical model is its ability to converge towards the exact 
solution as the number of mesh points is increased. For a source term m that is inde- 


pendent of time we can calculate the analytic steady state solution using Section 5.2.2 
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Experiment | Divide Thickness (m) | Boundary Position (km) 
Avg EISMINT 2856.9 594.9 
Moving Mesh 3005.3 576.9 
Exact 2951.5 579.81 


Table 5.3: Comparison between the average EISMINT experiment results, the moving mesh 


approach and the exact steady state solution. 


N_ | Boundary position (km) | Divide thickness (m) 
16 576.894 2995 
ol 579.104 2983 
61 579.636 2981 
Exact 579.814 2952 


Table 5.4: Table of 1D values for boundary position and divide thickness with varying number 


of mesh nodes. 


to provide exact values. We shall check convergence of the moving mesh scheme at the 
steady state using the EISMINT data (see Section 5.5.2), for an increasing number of 
mesh points. 

Table 5.4 shows that both the position of the moving boundary, Xy and the ice 
thickness at the divide, H; (= H(X,)) converge toward the exact solution as the num- 


ber of mesh points increases. 


Shown in Table 5.5 is the relative error of the boundary position and the relative 
error of the divide ice thickness with respect to the reference value. We see that the 
error decreases for increasing numbers of mesh points. Assuming that relative error 
for the boundary is x N~? where p is the order of convergence, we find p 2. The 
Lz norm of the relative ice thickness (discretely defined as H = (M1, Ho,...Hn)) shows 


evidence of superlinear convergence. 
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N 


| Xexact - X,| /| Xexact | 


| Hexact = H|o/|| HoX2" ||, 


16 
dl 
61 


0.0050 
0.0012 
0.0003 


0.0093 
0.0052 
0.0021 


Table 5.5: Error analysis for varying mesh nodes. 


We therefore find that the numerical CMF moving mesh method is convergent with 


second order accuracy. 
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Figure 5.6: Steady-state solutions to the EISMINT scenario: a) A 2D fixed grid method. 
b) 1D Moving mesh flowline radial model converted to 2D. c) Ice thickness profile along the 
flowline: exact (red), moving mesh (dash), fixed grid (solid). d) Diffusive velocity: moving 
mesh (dash), fixed grid (solid). 
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5.6 Summary 


In this chapter we began by looking at the shallow-ice glacier PDE and some of its 
features. In particular we looked at the diffusive component of velocity at the moving 
boundary. We developed a condition upon the profile of the ice thickness that deter- 
mines if there exists a finite diffusive velocity at the boundary or if that element is 
stationary. This condition was shown to be independent of temperature, density and 
the presence of a topographical bed. 

We then developed a procedure for modelling the flow using relative mass con- 
servation to determine a velocity comprised of the diffusive velocity and accumula- 
tion/ablation. This method was then approximated numerically by a moving mesh 
solution using finite differences, along with a modified method for a radially symmetric 
solution. 

The numerical model was then tested, firstly against a test function to demon- 
strate the model is able to achieve the advance, retreat and stationary movements that 
glaciers are capable of. Secondly we applied the model to the EISMINT moving mesh 
scenario, where the results compare favourably against both the fixed grid methods 


used previously and the exact solution. 


Chapter 6 


Data Assimilation on Moving 


Meshes in 1D 


In this chapter we look at applying the concepts of data assimilation to combine ob- 
servations with our moving mesh model to provide a best estimate of the real state of 
the ice sheet. We apply an optimised 3D-Var scheme to the model introduced in the 
previous chapter and show the impact on the components of the scheme when the mesh 
is moving. We test this scheme using a twin experiment applied to the test scenario in 
the previous chapter. 

We then proceed to build an extended scheme that treats the numerical mesh as a 
set of unknown state variables in order to include observations of positional features. 
We demonstrate how these observations alter our numerical predictions and highlight 


some of the challenges that arise with this extension. 


6.1 3D-VAR on a Moving Domain 


The observational network across the Earth has improved dramatically in recent years, 
giving us the opportunity to apply data assimilation to models in the cryospheric field. 
Here we apply the 3D-VAR scheme to our shallow ice moving mesh model. When 


assimilating on a moving mesh it is clear that the numerical domain will contain errors 
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in its representation of the physical domain, which does not occur in a fixed grid model. 
We now detail the differences in the algorithm when assimilating on a moving mesh. 

As described in Chapter 4 the optimal 3D-VAR analysis at time t, is given by the 
equation: 


Zi. = zi + Ki.(yx << C,2!), (6.1) 


subject to the dynamical model 


zi = Mazi._1, (6.2) 


where zz € R®@ is a vector containing the unknown state variables. The superscript ® 


represents the ‘analysis’ or best estimate of the state of the system and the previous 
guess or forecast is denoted using the superscript /. 
To obtain an analysis solution the forecast is corrected at each assimilation time by 
a weighted difference between the observations y; € IR? and the predicted observations 
C nz, where C, is known as the observation operator. The matrix weight K is optimally 
[65] chosen to be 
K, = BpC, (C, B,C; + Ry), (6.3) 


where B, € R?*? and R, € R®*? are the background and observation error covariance 
matrices respectively. 

The 3D-Var method is preferred to the other methods introduced in Chapter 4 as it 
provides an optimal explicit equation for the analysis solution. In addition the choice 
of the background error covariance matrix in the moving mesh environment can reduce 
the impact of some of the deficiencies of the 3D-VAR scheme (see Section 6.1.3). 

We now define the components of the 3D-Var scheme on the basis that the state 


vector contains the ice thickness variables only, i.e. 
Zy = Hy, (6.4) 


with z, € R%, for the N mesh points in the model. 
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6.1.1 Observation Operator 


The observation operator maps the state vector to observational space so that it is 
directly comparable to the observations. The primary variable in our shallow ice model 
is the thickness of the ice, h. Observing the ice thickness directly provides a value for 
the thickness itself, along with positional information for the location of the observation. 
Since the location is unlikely to align with the mesh points in a model, we use a linear 
interpolation technique to express the observation in terms of the nearest mesh points. 
An appropriate form for a single observation y; located between the mesh points 7;_1 


and x; would be 
yj = Pihist Pohi, (6.5) 


where 6; and (5 are interpolation coefficients defined as follows. Let observation y, 


occur at known location «*, then 


B, = —— ~~, and (6.6) 
Li Vi-1 
vu — vj 

es (6.7) 
Li Li-1 


Therefore for every observation of ice thickness, the corresponding row 7 in C contains 
the coefficients 8; and 85 located at Cj; and C;,; respectively. We assume that every 


observation lies within the numerical domain of the model. 


6.1.2 Observation Error Covariance Matrix 


The observation error covariance matrix R is a symmetric, positive definite p x p matrix. 
It represents the uncertainty in the errors in the observations which may arise from 
sources such as inaccuracies of the measuring instrument, subscale features that are too 
small to represent in the model or errors in the observation operator C. As this matrix 
is related to errors in the observations and not the dynamical model it is independent 
of the moving mesh method. 

A common assumption is to treat the errors in each observation as independent of 


the others, a reasonable assumption with direct measurements. Using observations of 
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the ice thickness means we may treat R as a diagonal matrix. If each observation has 


the same error variance, 07, we may write the observation error covariance matrix as 


R=o7I, ReR?*?. (6.8) 


6.1.3. Background Error Covariance Matrix 


The background error covariance matrix B is a symmetric, positive definite q x q matrix 
representing the uncertainty in the errors in the prior state of the system. 

This matrix is crucial, particularly in the optimal 3D-VAR scheme as it contributes 
towards the spread of information from observations to the state variables. It is common 
to approximate the error covariances using an analytic correlation function. There are 
many options available, the simplest would be to define B as a diagonal or tri-diagonal 
matrix; however, in a problem such as ice sheet modelling observations can be sparsely 
distributed and these simple options limit the spread of information from observations 
so are generally insufficient [28]. 

Instead we define a modified Gaussian function [91] which is defined over the entire 
domain and depends on the distance between the grid points of variables. The elements 
of B are given by 

big = opexp LOOM iF =1,..N (6.9) 


where L is the inverse of a background correlation length scale and o} is the error vari- 
ance associated with each ice thickness variable. On an evenly spaced grid each row 
(or column) forms similar functions across the domain centred on different values of x 
(see Figure 6.1(a)). 

The modified Gaussian function Eq. (6.9) is not restricted to an evenly spaced grid, 
which is beneficial to both h- and r- adaptive mesh methods. Moreover as the mesh 
evolves according to a moving mesh approach, the function reflects this with larger val- 
ues in locations where the mesh nodes have moved closer together and a reduction when 
they have moved further apart. For instance using the test scenario from Section 5.5.1 
with 26 mesh nodes and {3 = 1, when t = 12000 the same rows of the matrix now have 


the structure in Figure 6.1(b) since the nodes near the glacier front are closer together 
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Figure 6.1: A modified SOAR equation for rows of B relating to different places in the 


domain. a) evenly spaced grid with even structure, b) unevenly spaced mesh after evolution. 


than those near the divide. 

By viewing all of the rows together we are able to see the global change in the 
covariances as a result of the change in mesh. Figure 6.2 shows the elements of the 
background error covariance matrix, where the colour in position 7, 7 represents the size 
of B,,;. Initially the mesh is evenly spaced and the matrix B has a striped diagonal 
structure (Figure 6.2(a)). As the domain increases in size and the mesh nodes spread 
out, we see two main differences (Figure 6.2(b)). Firstly there are smaller covariances in 
the outer regions since the domain, and thus the overall distance between mesh points, 
has increased in size. Secondly, the mesh is no longer evenly spaced so the striped 
diagonal structure has curved to reflect areas where the nodes are closer, specifically 
towards the moving front. This allows control on how information is passed between 
the variables as the mesh evolves. In these figures the number of mesh nodes has been 
increased to 121 in order to demonstrate the changes smoothly. 


One of the drawbacks of the 3D-VAR method is that the background error co- 
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variances are static, whereas by defining the matrix B in the form above involves a 
recalculation at each assimilate time. The error covariances are then determined by the 
mesh evolution and implicitly depend on the non-linear flow of the numerical model. 
More advanced data assimilation schemes such as the Kalman Filter method intro- 
duced in Section 4.2.2 directly evolve the background error covariance matrix is using 


the model equations, generally in a linear form. 


1 
R 0.9 
B 0.8 
R 0.7 
Hl 0.6 
K 0.5 
R 0.4 
R 0.3 
100 
F 0.2 
: 120 0.1 
20 40 60 80 100 120 20 40 60 80 100 120 


(a) (b) 


Figure 6.2: Background error covariance matrix B using the Gaussian function Eq. (6.9); 
a) at the initial time when the mesh is evenly spaced, b) at t = 12000 when the mesh has 


moved. 


6.1.4 Analysis Solution and Mass Conservation 


By perfoming a step of an assimilation scheme the mass of the glacier is altered using 
external information in the form of observations of ice thickness. The CMF moving 
mesh method detailed in Section 5.3.1 is based on conserving relative mass. It is clear 


that the conservation principle (Eq. (5.42)) 


1 ee a 
xa | hoe = ia), (6.10) 


used to define the mesh deformation velocity and the ice thickness profile no longer 


holds when an assimilation step is made. Therefore an analysis solution or ‘reset’ of 
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some of the values is required before the numerical model can continue. Therefore at 
every assimilation cycle, after calculating the analysis solution to the ice thickness using 


Eq. (6.1), we recalculate the total mass @ such that the analysis solution of the mass is 


be (t) 
0“ (t) af he dae (6.11) 
0 
and the mass fraction constants jz are redefined as 
Tf de = un" (6.12) 
—— Se Sra), 6.12 
O°(t) Jo 


before the forecast step in the dynamical model. The analysis mass and equivalent mass 
fraction constants then correspond to the analysis estimation of the state variables. 

Alternatively we could treat the total mass variable 6 as an unknown variable and 
include it within the state vector (Eq. (6.15)). This would involve extending the back- 
ground error covariance matrix B and as a result the complexity, but would open up 
the possibility of including ‘observations’ of the total mass from sources such as energy 
balance or global climate models. 

The reset requirement is specific to the conservation based methods of mesh move- 
ment, models which do not depend on the total mass will not experience this issue. In 
this work we will calculate the analysis solution to the total mass and mass fraction 


constants at the end of the assimilation step. 


6.1.5 3D-Var Algorithm for CMF Moving Mesh 


We now present the 3D-VAR algorithm for assimilating on a moving domain using the 


conservation of mass fractions method from Section 5.3.1. 


1. Calculate a forecast of the state vector z/ by evolving the numerical model of 
the CMF method (Section 5.4.1) using the previous analysis solution as model 
inputs. This gives forecast values of the numerical mesh XS, the ice thickness H/ 


and the total mass @/. We also know the constant in time mass fractions ju. 


2. Use the optimal 3D-VAR formulation (Eq. (6.1)) to produce a new analysis solu- 


tion of the ice thickness H® using observations available at the current time. The 
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analysis solution to the numerical mesh is identical to the forecast, with X* = X/. 


3. Calculate the analysis solution of the total mass 6° and update the mass fractions 


wo" using Eqs. (6.11) and (6.12). 


4. Evolve the analysis solution using the dynamical model to produce a forecast at 


the next time step. 


5. Repeat steps 2-4. 


6.2 <A Test Experiment 


We will now test the 3D-VAR method applied to the moving mesh method using the 
test scenario in Section 5.5.1. Since we are using a test scenario, real observations do not 
exist and need to be generated from the numerical model. To achieve this we assume 
the true initial conditions are known and the model is perfect. The system is evolved 
with these initial conditions and the output recorded. Observations are then taken as 
a sample from this output, with random noise added to replicate observational error. 
We then choose a set of initial conditions that are incorrect, from which we evolve 
the system until a chosen time when we apply the data assimilation scheme using 


observations from the previous output. This procedure is known as a twin experiment. 
6.2.1 Initial Conditions 


Reference Solution 


To generate the observations for a twin experiment we first define the known set of 


exact initial conditions. Let the true initial ice thickness be the function 
h=1-27, xe (0,1). (6.13) 


The system is evolved using the method described in Section 5.4.1, with At = 0.2a and 
the initial domain spacing, Ar = 0.01m, corresponding to 101 mesh points. Both the 
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reference and forecast models evolve until an end time of T’ = 2000a, as we want to see 
the effects of including observations before the models reach a steady state solution. 


This design forms our reference solution, the output that provides our observations. 


Forecast Model 


For our experiments design we estimate the initial ice thickness by the perturbed func- 
tion 


h=1.1025—27, 2x € (0,1.05), (6.14) 


where we note that not only is the estimate of the ice thickness different but also the 
initial domain. As such the estimate of the total mass is different, along with the mass 
fractions used in the moving mesh method. The time step is taken to be the same as 
the reference solution so the observations can be taken at times that match the model 
without the need for temporal interpolation. We use 51 mesh points, an initial spacing 
of Ax = 0.021m. 

By using a higher resolution model for the reference solution we are able to generate 
more accurate values to use as observations. While this has the potential to introduce 


errors into the system, experiments have determined that these are negligible. 


6.2.2 Results 


We include sets of 19 observations of ice thickness from the reference solution, taken 
at t = 500 and t = 1500, sampled evenly across the domain with random noise taken 
from a normal distribution ¢, ~ N(0,02) with o? = 0.05m. The background error 
variances 0} is also set to 0.05m. This represents an error of approximately 5% of the 
corresponding variables. 

At the first assimilation time, t = 500, it can be seen in Figure 6.3(a) that there 
is a large difference between the forecast and reference thickness profiles. The analysis 
profile on the otherhand resembles the reference, except in the region at the front where 


the numerical model predicts there should be ice while in fact there is none there. This 
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is because there is no analysis solution to the numerical mesh. By the second assimila- 
tion time, t = 1500 (Figure 6.3(b)), we see the forecast profile, which has evolved from 
the analysis solution at t = 500, is closer to the reference profile. This suggests that 
by solely assimilating observations of ice thickness future evolutions of the dynamical 
model implicitly propagate information over subsequent time steps. After assimilating 
at this time we see the analysis profile is even closer to the reference solution. 

The implicit spread of information to the numerical mesh is emphasised in Fig- 
ure 6.3(c) where we see the moving boundary point x = b(t) being pulled towards the 
reference location after each assimilation time. 

Several experiments with a variety of initial conditions for x and h, along with vary- 
ing samples of observations produced similar levels of improvement to the estimate of 


the ice thickness profile. 


6.3 An Extended Scheme: Assimilating the Mesh 


Data assimilation methods have been devised for a variety of applications, with the 
general assumption that the physical domain from which observations are taken and 
the corresponding representation by the model are known exactly (or at least assumed 
to be known exactly). 

As we have seen the dynamical ice flow scenario (amongst others) occurs on a domain 
that is evolving in shape and size. In essence, the domain can be as much of an unknown 
as the variables within, especially since it is evolved using the dynamical model. With 
this in mind we propose an extension to existing data assimilation methods to directly 
include the unknown numerical mesh used in a moving mesh modelling process, such 
as the CMF method in Section 5.3.1. 

Of course, the mesh itself has no physical meaning as it is a numerical tool to solve 
the model equations. While the location of the majority of mesh nodes is arbitrary, there 
are some exceptions. In particular the domain boundaries are real, physical locations 
that are assigned mesh nodes. Similarly there could be other key moving features in the 


interior of the domain which could be assigned mesh nodes in a moving mesh scheme. 
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Figure 6.3: 1D moving mesh with data assimilation:- a-b) show the forecast, analysis and 
true ice thickness profile at t = 500 and 1500 respectively, indicating the observations after 


noise is added. c) shows the time evolution of the boundary point x = b(t). 


For example r-adaptive mesh methods are often used to locate the grounding line (see 
[39; 24]), the point at which grounded ice begins to float. If the locations features such 
as this are observable then we may wish to include these types of observations. In this 
work we focus solely on observing the boundary. 

In one-dimension the two mesh points at either end of the domain X, and Xy are 
situated on the boundary and represent a physical location. These locations, fixed 
or otherwise, are observable points that we may make use of in a data assimilation 


scheme. While there is a limited number of these ‘physical’ mesh nodes, most moving 
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mesh methods rely on a known connectivity between points and so the extension to the 
data assimilation scheme must include all the nodes within the domain to avoid mesh 


tangling. The state vector z in Eq. (6.1) becomes 


xX 
Fig = , (6.15) 
H 


to accomodate the mesh nodes. The use of subscript e determines a component that 
has been extended. 

The size of the state vector has increased by the number of mesh points (NV) in the 
model. Since our state vector previously only contained ice thickness variables located 
at each of these points, ze € R?% is now twice the size of the original state vector. 
Extending the state vector will alter the other components of the data assimilation 


scheme, which we shall detail here. 


6.3.1 Observation Operator Extension 


The observation operator for observations of ice thickness can be extended by writing 


Ce= (0 c) (6.16) 


where C,. € R?*?% and C € R®*% is the observation operator before the extension 
detailed in Section 6.1.1. 
In addition we wish to include observations of positional features, such as the bound- 


ary location. The observation operator may then be written 


C, 0 
C. = (6.17) 
aE 6 


where the rows in C, correspond to observations of positional features. These observa- 
tions differ from standard observations in that instead of each observation representing 
a variable with a corresponding position, the information is solely of the position of a 
feature itself. 

Since in one-dimension our moving mesh model estimates these features by asso- 


ciating a mesh point to them, there is a direct one-to-one relationship between the 
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observation of the position and the model estimation of this position. Therefore no in- 
terpolation is required and each row in C, consists solely of a 1 in the column associated 


with the mesh point for that boundary. 


6.3.2 Observation Error Covariance Matrix Extension 


The extended observation error covariance matrix still assumes the observations are 


independent, and as such we may write 


0 R 
Ree (6.18) 
o2,1 O 


where o2, is the error variance relating to observations of positional features. 


6.3.3. Background Error Covariance Matrix Extension 


The biggest change under the extended scheme lies within the background error covari- 


ance matrix. We may write the matrix in the form 


Bae ee ee (6.19) 
(B.,)' B 
where B,, represents the error covariance matrix between the mesh points and B is 
the previous matrix of error covariances for the ice thickness variables. The matrix 
B,, represents the covariances between errors in the mesh points and the ice thickness 
which we refer to as the cross-covariance matrix. 

The error covariances relating to the ice thickness, matrix B, are the same as the 
original background error covariance matrix defined in Section 6.1.3 using a modified 
Gaussian function (Eq. (6.9)). Similarly the errors in the mesh node locations can also 
be defined by their proximity to each other and as such we may define their covariances 


using the same modified Gaussian function 


{Box} ij 7 Cnr exp P@@O-20))" ce) = 1, el (6.20) 
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with a different variance op. 

Since the mesh is an artificial representation of the domain, it is not immediately 
clear how to define B,;,. While there exists a relationship between the mesh nodes and 
the ice thickness variables through the moving mesh conservation principle (Eq. (5.42)), 
it is non-linear. Other moving mesh methods may contain a linear relationship, but for 
the method used in this thesis we seek an alternative. One approach is to assume there 
is no relationship at all and define the cross-covariance error matrix as a N x N zero 
matrix 


B,, = 0. (6.21) 


This means that at the point of assimilation no information from observations regarding 
the mesh is passed to the ice thickness variables and vice versa, with the matrix B 
taking the form in Figure 6.4(a). Therefore unless there are observations of both the 
ice thickness and the mesh nodes, the analysis solution relies on future progression of the 
numerical model to filter information through and the extended scheme is superfluous. 

A second approach relates the mesh nodes that correspond to physical points with 
the other variables. In our model the first and last mesh point correspond to the 
physical domain boundaries, which we relate to the ice thickness variables across the 
domain. We define this relationship using the modified Gaussian function Eq. (6.9), 


where the corresponding rows of B,,;, are given by: 


{Ban} 1; = hep exp HOO", and (6.22) 


Bint y; = % exp L@N)-2G)” a emer s 6.23 
Nj bah 


with variance o7,,. 

The structure of the matrix B under these two approaches is visible in Figure 6.4 
and we shall discuss the impact of each method through the experiments below. 

A major drawback of extending the scheme is the size of the background error 
covariance matrix, which is four times the size. For systems with a large number of 


mesh points this may increase the computational requirements to an unfeasible level. 
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Figure 6.4: Background error covariance matrix for the extended scheme; a) without cross 
covariance terms between the mesh and the ice thickness, b) with cross covariances between 


the boundary points and the ice thickness providing additional rows/columns. 


6.3.4 Extended 3D-Var Algorithm for CMF Moving Mesh 


The algorith for the extended 3D-VAR scheme applied to the CMF moving mesh method 


is given by: 


1. Calculate a forecast of the state vector z/ by evolving the numerical model of 
the CMF method (Section 5.4.1) using the previous analysis solution as model 
inputs. This gives forecast values of the numerical mesh XS, the ice thickness H/ 


and the total mass @/. We also know the constant in time mass fractions ju. 


2. Use the optimal 3D-VAR formulation (Eq. (6.1)) to produce a new analysis so- 
lution of both the numerical mesh X° and ice thickness H® using observations 


available at the current time. 


3. Calculate the analysis solution of the total mass 6° and update the mass fractions 


pe” using Eqs. (6.11) and (6.12). 


4. Evolve the analysis solution using the dynamical model to produce a forecast at 


the next time step. 


5. Repeat steps 2-4. 
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6.4 Extended Scheme Experiments 


Now that we have an extended data assimilation scheme capable of producing an anal- 
ysis solution to the numerical mesh we will use the experimental design in Section 6.2 
to investigate this scheme and its ability provide a better estimate of the state of the 
system. The observation noise is kept constant between these experiments so results are 
directly comparable. We will evaluate the impact of using observations of the bound- 
ary within the scheme, both with and without the cross covariance terms introduced in 
Section 6.3.3. The results are compared to the solution before the extension, as seen in 


Figure 6.3, which we shall refer to as the benchmark solution. 


6.4.1 Observing the Boundary 


We now include an observation of the moving boundary (x = b(t)), in addition to the ice 
thickness observations, at both assimilation times, t = 500 and t = 1500. The boundary 
observations are taken from the reference solution with random noise sampled from a 
normal distribution €, ~ N(0,02,) with o2, = 0.05m. To demonstate the impact of in- 
cluding these observations let the cross covariance background error variance 07,,,, = 0m 
(i.e. Eq. (6.21) holds), which means no information is passed from the thickness profile 
to the mesh at the time of assimilation. As you can see in Figure 6.5(a), the extent of 
the domain is corrected along with the thickness profile. Since no information is passed 
between the mesh and the ice thickness, at the time of assimilation the analysis profile 
is actually a transformed version of Figure 6.3(a), i.e. the profile is identical but on a 
squeezed domain. 

In Figure 6.5(b) we notice that the forecast profile is significantly closer to the true 
ice thickness profile than the forecast profile without a boundary observation shown in 
Figure 6.3(b). We can also see a ‘jump’ towards the true location of the boundary at 
the time of assimilation in Figure 6.5(c). As a result the convergence towards the true 
boundary location is significantly quicker than before. 


It is also vital that the mesh points do not tangle when we alter the mesh at the as- 


CHAPTER 6. 1D DATA ASSIMILATION 90 


—— Forecast 

— Analysis 
—Truth 

+ Ice Thickness Obs 
* Boundary Obs 


Ice Thickness (m) 


(a) 
15 
ee - . ——Forecast 
i eS — Analysis 
i ee —Truth 
L i a * Ice Thickness Obs 
+ Boundary Obs 


Ice Thickness (m) 


06 08 
Distance from Divide (m) 


145 ae 
1b 
Ey ae 
Ss ——— 
5 val — —— 
B 125 as ee 
Qa | = 
Sa = = 
GB yisL eee 
2 —— Forecast 
3 1 bh — Analysis 
1.05 /- SS —Truth 
‘oe 
re ! | ! | | | | | | | 
r) 200 400 600 200 7000 1200 7400 7600 7800 2000 
Time (a) 


Figure 6.5: 1D moving mesh with data assimilation including observations of the boundary:- 
a-b) show the forecast, analysis and true ice thickness profile at t = 500 and 1500 respectively, 
along with the observations after noise is added. c) shows the time evolution of the boundary 


point x = b(t). 


similation point. Figure 6.6 shows the evolution of the mesh points for the benchmark 
scenario and the case including observations of the boundary. Naturally the bench- 
mark solution evolves regularly with no influence on the mesh from any observations 
(Figure 6.6(a)), while the impact including a boundary observation in the assimilation 
is visible in Figure 6.6(b) by a change in direction of mesh nodes near the boundary. 
When the boundary observation is included we see an immediate knock on effect on 


the interior mesh points all the way across the domain, evident by the sharp ‘ridge’ 
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Figure 6.6: 1D moving mesh with data assimilation:- evolution of mesh nodes through time. 


a) benchmark solution, b) including observations of the boundary. 


that appears at the assimilation times. This is due to the Gaussian function in B,, 


(Eq. (6.20)) spreading the information of the boundary location through the rest of the 


mesh. 


6.4.2 Including Cross-Covariances 


This next experiment analyses the impact of the cross-covariance part (Eqs. (6.22) 
and (6.23)) of the background error covariance matrix B detailed in Section 6.3.3. We 
set a}, = 0.0075m in Eqs. (6.22) and (6.23) to allow information to be passed between 
the mesh and the ice thickness, but with a weaker relationship than exists between the 
mesh nodes or the ice thickness variables. 

We see in Figure 6.7(a) that without observing the boundary the analysis estimate 
of the domain is an improvement over the forecast. While the improvement is not 
as significant as the scenario with a boundary observation (Figure 6.5(a)), it is still 
preferable to the benchmark solution where the domain is fixed at an assimilation time. 

By including both the cross-covariance terms and a boundary observation we are 
making use of all the information available, which results in a further improvement at 


the first assimilation time, t = 500, as can be seen in Figure 6.8(a). By the second 
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Figure 6.7: 1D moving mesh with data assimilation including cross-covariance terms in the 
extended background error covariance matriz, Be:- a-b) show the forecast, analysis and true 
ice thickness profile at t = 500 and 1500 respectively, along with the observations after noise 


is added. c) shows the time evolution of the boundary point x = b(t). 


assimilation time, t = 1500, the forecasted profile is virtually indistinguishable from 
the true profile. At this point the random noise on the observations means the analysis 


solution is a slightly worse estimate of the truth, though the amount is negligible. 


6.4.3 Discussion 


We have shown that by extending the data assimilation scheme to incorporate the 


numerical mesh the prediction of the numerical domain and as a result the overall ac- 
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Figure 6.8: 1D moving mesh with data assimilation including cross-covariance terms in B 
and boundary observations:- a-b) show the forecast, analysis and true ice thickness profile at 
t = 500 and 1500 respectively, along with the observations after noise is added. c) shows the 


time evolution of the boundary point x = b(t). 


curacy is improved. The largest improvements in accuracy arose when observations of 
the boundary were directly included, while defining a cross-covariance matrix to spread 
the information from observations between the mesh and the ice thickness also saw 
improvements over the ordinary scheme. 

In addition, experiments designed to explore the effects of the quantity of observa- 
tions, along with their distribution in space and time show the scheme, both ordinary 


and extended, consistently improve the overall accuracy and are robust enough to cope 
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with the variations. 
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6.5 Summary 


In this chapter we presented the 3D-VAR data assimilation scheme and detailed the 
differences to the various components that make up the scheme when applied is a mov- 
ing mesh compared with a fixed grid. We discussed the need for an analysis solution to 
the total mass and mass fractions for the chosen method of mesh movement to ensure 
the mass conservation required. 

We then tested the scheme using a twin experiment, where results showed the scheme 
can achieve improved predictions of the ice thickness profile at the moment of assim- 
ilation. This improvement then implicitly improved the estimate of the domain over 
future evolutions of the model. 

We then provided an extension to the scheme to include the numerical domain in 
the state vector, with the aim of improving our prediction of the domain at the as- 
similation time. This allowed observations of positional co-ordinates to be included. 
The extended scheme included changes to the structure of many of the features in the 
3D-VAR methods which we detailed. 

We saw the benefits of the extended scheme to the analysis solution. We examined 
the impact of including an observation of the boundary, along with a cross-covariance 
relationship between the errors in the numerical mesh and the ice thickness profile. 
By including these cross-covariance terms we showed the improvement to the domain 


prediction, both with and without observations of the boundary. 


Chapter 7 


A 2D Moving Glacier Model 


In this chapter we look at the shallow ice model in two horizontal dimensions. We 
begin by considering the total mass of the glacier and how it evolves over time, before 
non-dimensionalising the shallow ice PDE to remove the physical units. 

We then extend the CMF moving mesh procedure for modelling ice sheet flow into two 
spatial dimensions by deriving a weak formulation of the relative mass conservation 
principle. The weak form may then be approximated using a finite element approach. 
We demonstrate the method using a test problem, before comparing results of the 


EISMINT scenario to both fixed grid and exact solutions. 


7.1 The 2D Shallow Ice PDE 


A one-dimensional flowline model such as the one presented in Section 5.3 is useful for 
providing a general view of simple glacial movement, but it is only applicable in certain 
scenarios. We therefore extend this model into two dimensions to include the flow in 
both horizontal directions with the aim of producing a more complete picture of the 
movement. In two dimensions the shallow ice mass balance equation (Eq. (2.14)) takes 


the form 


hy =m — V.(hu) (7.1) 
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Figure 7.1: Arbitrary two-dimensional domain Q(t), with a moving boundary exhibiting a 
Dirichlet condition at the glacier front (Up) and Neumann boundary conditions (Un) along 


the rock walls 


where h = h(x,t) is the thickness of the ice and m = m(x,t) represents the ice- 
equivalent accumulation rate, both given in terms of positional coordinates (x,y). The 
diffusive velocity u is defined to have components (u’,u”)’, which describes the flow 
in the x and y directions respectively. 

A two-dimensional domain, Figure 7.1 for example, possesses two types of boundary. 
In the upper regions of a glacier and along the valley sides the boundary is defined by 
a rock face or slower moving ice, which we denote by 'y. There can be no flow across 
this boundary, so there exists is a no-flux Neumann condition written as 


Ooh 

an 0, un=0 onTy. (7.2) 
In the regions not bounded by these geographical features the ice meets the glacier bed 
and the boundary is free to move. This (Dirichlet) boundary is denoted [p with the 


condition 


f=0 onl pg: (7.3) 
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7.1.1 Integral Property of Total Mass 


We first consider the total mass of the glacier and its evolution over time. As in the 
one-dimensional case let @ represent the total mass of the system. In two dimensions 


on a domain 2(t) this is defined as 


9 = i hdx, (7.4) 
(8) 


under the assumption that density is constant. Differentiating Eq. (7.4) with respect 
to time using the Reynolds Transport Theorem (see Section B.4 for details, [90]) gives 


es hdx, (7.5) 
dt Q(t) 
= . hydx + f (van)hds. (7.6) 
Q(t) E 
where v = |v(x,t)| is the component of the velocity at the boundary. Under the 


boundary conditions Eqs. (7.2) and (7.3) the surface integral term is zero. This leaves 
a volume integral over the whole domain. Substituting for h; from the mass continuity 


PDE, Eq. (7.1), gives 


j= mdx— f V.(hu)dx. (7.7) 
Q(t) Q(t) 


In order to evaluate the divergence in the second term of Eq. (7.7) we note that since 
(hu) is continuous and differentiable we may apply the Divergence Theorem (see Sec- 


tion B.5) to obtain 
6 = mdx — f (i) nds (7.8) 
Q(t) ij 
where again application of the boundary conditions force the surface integral to be zero. 


This leaves us with the expression 


at.» 
6= — hdx = mdx (7.9) 
dt Jor) a(t) 
to represent the change in total mass. As with the one-dimensional case this means 
that any change in the ice mass comes from the net effect of the source term over the 


entire glacier. 
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7.1.2 Non-Dimensionalisation of the Shallow Ice PDE 


For many problems it is beneficial to non-dimensionalise the variables in order to avoid 
extreme values and understand the relative strengths of the features. Each variable 
within the PDE is scaled by a constant parameter to leave a dimensionless variable. 
The transformed variables for ice thickness, accumulation and flow velocity are given 
by 

(7.10) 


where the non-dimensional variables are represented with an overbar and the scaling 


parameters by square brackets. Similarly the non-dimensional time and space variables 


are 
- ¢ x 
t=—,x=—. (ld) 
?] !] 
We can therefore write the shallow ice PDE from Eq. (7.1) in the form 
ae = |[m|m — lv. (7.12) 


It is important that the scaled shallow ice expression, Eq. (7.12), maintains the same 


balance of terms as the original PDE. Therefore to maintain balance 


(A) ay — ele 
Gi o 


By definition of the shallow ice equation, the flat bed flow velocity u = ch*|Vh|?Vh 


(7.13) 


is given in terms of the ice thickness and gradient, which implies the scaling for the 


diffusive velocity is of the form 


_ [Al 
[tel ve (7.14) 


Putting this into Eq. (7.13) implies that balance will be maintained provided 
inl" (e) = (a (7.15) 


Therefore to non-dimensionalise the shallow ice equation the user must define two 
scaling parameters from [h], [tf] and [l|. Then the remaining parameter in Eq. (7.15) 


and |m] in Eq. (7.13) can be found such that the non-dimensional shallow ice PDE in 
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Eq. (7.12) has the same balance of terms as the full equation. We can then solve the 
modified PDE using the same techniques as before, with the output re-dimensionalised 
at the end. To save notation becoming cumbersome we shall express equations hereafter 
using the original notation; however, the results shown were calculated using a non- 


dimensionalised model. 


7.2. Glacier Movement 


In one dimension we developed a moving mesh method to describe the net glacier move- 
ment using a conservation principle to maintain relative mass over time. This method 
was then implemented numerically using a finite difference scheme to approximate the 
relevent equations. For a two-dimensional moving problem, finite differences are not a 
viable option as it is difficult to estimate the diffusion terms in Eq. (7.1) and the risk of 
mesh tangling is high. Instead we seek a finite element approximation which admits a 
locally once differentiable function. With this in mind the CMF method can be derived 


in a weak form to facilitate the finite element approach. 


7.2.1 Conserving Mass Fractions (CMF) - Weak Formulation 


The relative mass of a glacial subdomain with constant density is conserved through 


time according to the CMF principle 


1 
al hdx = p, (7.16) 
0 w(t) 


where jz is constant in time and is applied to any moving subdomain w(t) € Q(t). The 
total mass @ is defined in Eq. (7.4). To express Eq. (7.16) in a weak form we replace 
w(t) by Q(t) and include a continuous, differentiable test function w, 


;[ whdx = pL, (7.17) 
O Jor) 
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where the w form a partition of unity. With this weak form we extract the net velocity 


by first differentiating Eq. (7.17) with respect to time to give 


sal (| wha) = yO. (7.18) 
dt Q(t) 


To perform the time differentiation on the left hand side of Eq. (7.18) we employ the 
Reynolds Transport Theorem (Section B.4) to change from a Lagrangian perspective 


to an Eulerian one to produce 


i (wh)dx+ (wh)(v.n)dP = p08, (7.19) 
Q(t) 


T(t) 


which by the product rule of differentiation gives 


| wrh + whydx + ¢ (wh)(v.n)d0 = p08. (7.20) 
Q(t) 


T(t) 
The boundary integral contains the normal component of the net velocity v that we 
seek on the surface elements dl’. In order to calculate the velocity v we convert the 


surface integral to a volume integral using the Divergence Theorem, giving 
i. we + wh; + V.(whv)|dx = p06, (i210) 
t 
where again utilising the product rule of differentiation 
is A (w, + v.Vw) + wh, + wV.(Av)|dx = p06. (7:22) 
t 


By making the assumption that w is frozen in the moving domain ()(t), which moves 


with velocity v, it follows that the test function must satisfy the advection equation 
w,zt+tv.Vw = 0. (7.23) 
Combining Eqs. (7.22) and (7.23) leaves 


| [why + wV.(hv)|dx = p06. (7.24) 
Q(t) 


We can substitute the rate of change in ice thickness h; from the shallow ice PDE, 


Eq. (7.1), into the first term of Eq. (7.24) to give 


| w(m — V.(hu)) + wV.(hv)dx = 8, (7.25) 
Q(t) 
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leading to 


: w¥. (hvac = 1d — | wma + | wV.(hu))dx. (7.26) 
Q(t) Q(t) Q(t) 


The net velocity v that we seek is now expressed in Eq. (7.26) in terms of spatial 
gradients. To extract the velocity itself we use integration by parts. Evaluating the left 


hand side of Eq. (7.26) gives 
| wV(vjdx= f (whv) a Vw.(hv)dx, (7.27) 
Q(t) 
aa Vw.(hv)d (7.28) 


after applying the boundary conditions, Eqs. (7.2) and (7.3), which force the boundary 
integral to be zero since either h = 0 or v.n = 0 on the boundary. 
Similarly the last term on the RHS of Eq. (7.26) may also be evaluated using integration 
by parts: 
i wV.(hu)dx = ¢ (whu).ndx — Vw.(hu)dx, (7.29) 
Q(t) Q(t) 


= Fl Vw.(hu)dx. (7.30) 
a(t) 


Again, enforcing the boundary conditions reduces the boundary integral in Eq. (7.29) 
to zero. Putting Eqs. (7.28) and (7.30) back into the expression for the net velocity 
(Eq. (7.26)) gives 

Vw.(hv)dx = | 


wmdx — p60 + i Vw.(hu)dx. (7.31) 
Q(t) 


Q(t) 
Eq. (7.31) is a weak formulation satisfied by the net velocity v as a combination of the 
movement induced by both the source and diffusion parts of Eq. (7.1) for an arbitrary 
point in a region. Any solution to Eq. (7.31) determines the mechanism by which the 
domain evolves over time, which we may approximate using finite elements. This can 
also be simplified to the one-dimensional case where results are comparable to the finite 
difference approach in Section 5.4.1, but that is not done here. 
We note that there is no unique solution to this weak formulation, since it is possible 


to add to Av an arbitrary velocity V x q for some vector field q in Eq. (7.26). Here we 
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assume that v is irrotational (not unreasonable for glaciers) and write the velocity in 
terms of a velocity potential w: 


ve Vy, (7.32) 


which is unique to a constant. Eq. (7.31) can then be written as 


if Vw.(hV)dx = i wmdx — pO + Vw.(hu)dx. (7.33) 
Q(t) Q(t) 


Q(t) 


7.3 A Finite Element Scheme 


Now that we have a weak form for the net velocity (in the form of a velocity potential, 
Eq. (7.33)), based upon the relative mass conservation principle we seek a finite element 
numerical approximation to model the shallow ice equation. 

Any domain (Q(t) can be discretised into a set of elements defined by a set of nodes 
X1, Xo...Xy where X; = (X,Y);. These nodes can be connected to form a mesh made 
of triangular elements across the domain. Let us define the test function w as a piecewise 
linear basis function, 


w = (x, t), (7.34) 


which takes the value 1 at node X; and zero at all other nodes, thus forming the 
‘pyramid’ function shown in Figure 7.2. The set of functions, @; has the property that 


they form a partition of unity 
N 
S d=) (7.35) 
i=1 


To ensure this partition when a Dirichlet boundary condition is imposed, the mass is 
shifted to the interior nodes. This means the pyramid function is altered for elements 
connected to a Dirichlet boundary (for more details see [47]). 

Choosing the basis functions in this way enables us to define a piecewise linear approx- 


imation of h(x) such that 


N 
he H= 5 ~ Hyd; (7.36) 
j=l 
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Figure 7.2: Two-dimensional linear finite element ‘pyramid’ basis function [109] 


where H, is the approximation to h at X;. Similarly M, is the discrete approximation 
to m(X,) with 


N 
ma M= 5° Mj¢;. (7.37) 


j=l 
The total mass @ (Eq. (7.4)) can be approximated by substituting H for h, then 
summing over all the triangles which make up the domain. An approximation © to @ 


is then given using the two-dimensional equivalent of the trapezium rule, 


P 3 
1 
oi Ss" Ap (>: Ha) (7.38) 
p=1 g=1 
where P is the number of triangles in the domain ()(t) and A, is the area of triangle p. 


Hpq is the approximation to the ice thickness at vertex q of triangle p. 


Similarly we may approximate the differential of the total mass, 6 (Eq. (7.9)), by 
iv 3 
C= S "Ap (>: My) (7.39) 
p=1 q=l 
where M,, is the approximation to the ice equivalent accumulation rate at vertex q of 


triangle p. 


With these approximations we can formulate a numerical expression to the relative 


mass fraction constants ju; (Eq. (7.17)). Using the basis function Eq. (7.34) and the 
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approximation of h from Eq. (7.36), the expression for can be written 


= ae whdx (7.40) 
all 

= 6.0 H, (7.41) 
s= 

= MH (ine) 


where we introduce the matrix M, the standard mass matrix for linear basis functions 


used in finite element methods with elements 


My = f didjAdx. (7.43) 


Using a similar approximation to that used for the ice thickness and the source term in 
Eqs. (7.36) and (7.37), respectively, we may approximate the velocity potential « by 
N 
p= So bibs, (7.44) 
j=l 
where w; = w(X;). The gradient of the velocity potential is therefore approximated by 
N 
Vb => ¥;Vo;. (7.45) 
j=l 
The initial ice thickness, accumulation term, total mass, mass fraction constants and 
velocity potential have now been approximated, which allow us to calculate the numer- 
ical estimations in Eq. (7.33) and thus find the net velocity. The algorithm for a single 


time step can be split into three stages. 


Stage 1 
In the first stage we aim to approximate the net velocity, for which we require the 
velocity potential from Eq. (7.33). 

Substituting the approximations Eqs. (7.32) and (7.45) into the equation for velocity 
potential, the left hand side of Eq. (7.33) can be written as 


: : Vw.(hVap)dx dX i ; 11,6. ox w; (7.46) 


= Kya (7.47) 
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where w is a vector with entries ~;. Kp is the stiffness matrix used in finite element 
methods weighted by the ice thickness h. The elements of the weighted stiffness matrix 
are given by 


Ap 


The velocity potential ~ is then found by a solution of the approximation to Eq. (7.33) 
) 
Knw = MM — QMH +f (7.49) 


using the approximations to h, m and pw from Eqs. (7.36), (7.37) and (7.42), with 
fi = J ¢:.(hu)dx. 


Finally, to recover the net velocity v we use an approximation to the weak form of 


Eqs. (7.32) and (7.45) to obtain 


N N 
es n ods Vi pe [fe enV Wj, (7.50) 
N 


j=l j=l 


MV 


2» . enV Yj. (7.51) 


j=l 
Stage 2 
In the second stage we use the calculated velocity V to move the nodes of the mesh by 
the explicit Euler scheme: 

Xe KP AV (7.52) 
where k denotes the time discretisation level. Similarly, the total mass O is also updated 


by an explicit Euler scheme using the approximated value of 6 from Eq. (7.39), 
of! — EF + Ato. (7.53) 


The time step is again based upon the general rule for diffusion equations, as in 
Eq. (5.63). 

Stage 3 

The final step is to recover the ice thickness at the new nodes. Rearranging Eq. (7.42) 
gives 


HH! = [Mtoe] a (7.54) 
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determines the ice thickness at the next time step while enforcing relative mass conser- 


vation. 


7.4 Experiments 


We now demonstrate the two-dimensional moving mesh method detailed in Section 7.3. 
As with the one-dimensional version of the method we begin with a test scenario to 
demonstrate the method, before again using the radially symmetric EISMINT scenario 
to compare the method with existing methods. Finally we will demonstrate domains 


that are not radially symmetric. 


7.4.1 Testing the Model 


For an initial test we shall use a circular domain centered at the origin (0,0) with radius 


b(t) and an initial profile of the form: 
hy = (6(0)? — |x|?)*. (7.55) 


The parameter a influences the profile and the diffusive velocity in a similar way to 
the analysis in Section 5.2.1. Initially the radius of the domain is given as b(0) = 1, 
where the boundary is entirely Dirichlet with the zero ice thickness condition given in 
Eq. (7.3). This means that the entire boundary is free to move. We choose a source 


term that is radially symmetric, i.e. 


a pss 


where 6 determines the equilibrium line and y is a parameter to control the scale. Using 
a radially symmetric domain and source term forces the model to behave evenly in all 
directions, which is ideal for testing the model. The source function is independent of 
time which means there exists a steady state solution. This may be found by using the 
radially adjusted version of the steady state solution, found in Eq. (5.82). It can be 


shown therefore that the steady state radius is given by 


bss = +28. (7.57) 
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Using the parameters defined in Table 7.1 the two-dimensional numerical moving mesh 
model is tested on a Strang type mesh (see Section A.2) with 177 mesh nodes, as shown 
in Figure 7.3. We will now demostrate the ability of the two-dimensional numerical 


model to simulate advance, retreat and an initially stationary front. 


Physical Parameters 


es Flow-law exponent 

A= 107%(Pa)%a-* Flow-law parameter 
g=9.81ms~? Acceleration due to gravity 
p = 910kgm-? Ice density 

c= —2Ag"p"/(n + 2) Constant Parameter 

y = 0.0005 Scale of accumulation rate 


Computational Data 


N = UGT Number of gridpoints 
At = 26 Time Step 
T Final Time 


Table 7.1: Values used within the test scenario including physical parameters of the PDE 


and the data used in the computational domain. 


Glacier in Advance: a = 3/7, 6 = 0.72, T = 30000a 

We first set the equilibrium line parameter to 6 = 0.72, so that the steady state 
solution has a larger extent than the initial domain. We see in Figure 7.4 that the ice 
thickness builds up over time across the domain, primarily close to the centre where the 
accumulation term in Eq. (7.56) is greatest. The boundary of the domain increases sym- 
metrically over time until the domain reaches steady state at bs, = V/1.44 in Eq. (7.57). 
Note that in steady state there still exists a small amount of movement in the interior, 
but the effect on the ice thickness over the domain is negligible and crucially there is 


zero velocity on the boundary which means the domain is stationary. 
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Figure 7.3: Initial nodal location and mesh, using the Strang method (see Section A.2) for 


the two-dimensional test scenario 


Glacier in Retreat: a = 3/7, 6 = 0.32, T = 10000a 

Now let us define the equilibrium line by 6 = 0.32, such that the extent of the steady 
state solution in Eq. (7.57) is inside the original domain. We observe in Figure 7.5 that 
the ice thickness still builds up in the center of the domain where the accumulation 
term is positive, only now the boundary retreats back to b,, = 0.64. Due to the sharp 
gradient of the accumulation term (Eq. (7.56)) the simulation reaches the steady state 
solution quickly as this is the dominant term in the net velocity Eq. (7.31). Encour- 
agingly, the numerical mesh nodes respond evenly in the domain, maintaining their 


connectivity and avoiding mesh tangling. 


Stationary Front: a=1, 6 =1, T = 2000a 

A final test of the model sets the equilibrium line by choosing 3 = 1 so that the accu- 
mulation term is zero on the initial boundary. In addition, the initial profile Eq. (7.56) 
is set with a = 1 so the analysis in Section 5.2.1 is applicable and the flow velocity is 
initially zero on the boundary. 

We see in Figure 7.6 that the net velocity remains zero around the boundary as a 
result of the choice of initial conditions, until the ice thickness profile builds up to meet 


the condition in Section 5.2.1, and the flow velocity can push the boundary forwards. 
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Figure 7.4: Advancing Glacier: Ice Thickness Profile (left) and Net Velocity (right) at 
t = 0, 2000, 15000, 30000 (in descending order) with a = 3/7 and 8 = 0.72. 
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Figure 7.5: Retreating Glacier: Ice Thickness Profile (left) and Net Velocity (right) at 
t = 0, 1000, 5000, 10000 (in descending order) with a = 3/7 and 8 = 0.32. 
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Figure 7.6: Initially Stationary Glacier: Ice Thickness Profile (left) and Net Velocity (right) 
at t = 0, 1000, 2000 (in descending order) with a= 1 and B=1. 
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7.4.2 European Ice Sheet Modelling INiTiative 


We now compare the two-dimensional moving mesh model with the EISMINT equations 
encountered in Section 5.5.2. As before we replicate the radially symmetric moving 
margin experiment and compare the results with the fixed grid methods presented in 
[49]. 


The source term is given as a function of the distance from the origin, 
m = min{0.5, 7(8 — |x|)}, (7.58) 


with 6 = 4.5 x 10°km determining the equilibrium and y = 10-2ma~'km7! the scal- 
ing factor. The experiment starts with zero ice thickness, and is evolved forward by 
one fixed grid Euler time step to get an initial profile for the ice thickness, whose dis- 
tributed mass we will preserve under the two-dimensional weak form CMF method in 


Section 7.2.1. Therefore the initial ice thickness profile is 
h® = At x min{0.5, y(8 — |x|)}, (7.59) 


over a circular domain with radius 6. The initial radius corresponds to the equilibrium 
line as ice has accumulated in this region. The steady state boundary is the same as in 


one dimension (Eq. (5.83)) which applies to the whole circumference of the domain: 
b., = 579.81km. (7.60) 


The physical parameters used in the experiment are given in Table 7.2 along with 
the numerical data. The initial mesh discretisation uses the Strang approach (see 
Section A.2) with 721 mesh nodes. The time step is taken to be 2a, with the simulation 
running to a final time T’ = 30000a. 


To assess the results of the moving mesh approach we directly compare the output 
with that of a fixed grid scenario discretised evenly on a 31 x 31 square, a total of 961 


grid points. We see in Figure 7.7 that the moving mesh approach resembles the circular 
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Value Quantity 

C3 Flow-law exponent 
A=10-" Pa ea-* Flow-law parameter 

g=9.81ms~? Acceleration due to gravity 

p = 910kgm-? Ice density 

m = min{0.5, (8 — |x|) } Ice equivalent accumulation/ablation 


=k Slope of source function 


y = 10-2ma-!km 
B=4.5 x 10°km Distance at which source function becomes neg- 


ative 


Table 7.2: Two-dimensional EISMINT model parameters and functions. 


nature of the domain significantly better than the fixed grid method despite having 25% 
fewer grid points. The boundary is not perfectly circular due to the mesh representation, 
though the jagged edges are less prominent than in the fixed grid case. More specifically, 
the mesh points that do lie on the boundary provide a direct approximation to the 
location, without the need for interpolation. 

Comparing both methods to the exact solution, we see in Table 7.3 that the moving 
mesh approach approximates the solution to the divide thickness to within 1%, while 
the error in the fixed grid method is twice as much. In particular the moving mesh 
adept at estimating the boundary position which it approximates to 0.14% of the exact 
location, whereas the fixed grid method has an error of 2.6%. When considering the 
scale is in kilometers this value is significantly higher. 

However, the moving mesh method took approximately 4 times as long to run on a 
standard desktop than the fixed grid method, despite having fewer grid points. This 
could probably be reduced using more efficient coding techniques, though some speed 


is invariably lost due to the overheads in finite element methods. 


Choice of Initial Mesh 
The initial mesh can potentially be important in achieving high levels of accuracy. 


To that end we introduce four different methods for generating a circular mesh, detailed 


CHAPTER 7. A 2D MOVING GLACIER MODEL 115 


2500 


1000 


500 


Figure 7.7: Steady state solution to the EISMINT scenario: a) A 2D fixed grid method, b) 


2D Moving mesh model. 


Fixed Grid | Moving Mesh 
Divide Thickness Error 1.96% 0.86% 
Boundary Position Error 2.60% 0.14% 
Computational Time 445s 2057s 


Table 7.3: Error comparison of steady state results from Fired Grid and Moving Mesh 
solutions to the 2D EISMINT problem. 


in Appendix A, to use for comparison. By comparing different meshes we may identify 
common themes of the underlying method. We will apply the EISMINT scenario to 
each of these initial meshes, with approximately the same number of mesh nodes. 

As expected the four meshes give comparable results (Figure 7.8) with only minor 
differences. The results from the hex mesh (Figure 7.8(c)) looks the most different from 
the others as this mesh has a large spatial gap between the boundary nodes and their 
connections to the interior, which results in long, thin triangular elements that skew the 
shading of the plot. As a consequence of having fewer mesh points near the boundary 
this mesh is the least accurate in resolving the boundary distance, as seen in Table 7.4. 


A common theme for all the meshes is to overestimate the thickness at the divide. 


CHAPTER 7. A 2D MOVING GLACIER MODEL 116 


This suggests that the flow velocity near the centre is slightly underestimated, as this 


acts to move ice towards the boundary. 


x 10° 


(a) Spoke 


x x10° 


(c) Hex (d) GMSH 


Figure 7.8: Steady state ice thickness for the EISMINT problem using four different initial 


meshes. 


Convergence 

To test the performance of the numerical models we increase the resolution and as- 
sess convergence towards the exact solution. We check convergence for the spoke mesh 
(Appendix A.1) using the EISMINT data found in Table 7.2, while alternate initial 


meshes and inputs exhibited similar behaviour. As you can see in Figure 7.9, as the 
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Mesh N_ Boundary (km) Divide Thickness (m) 
Spoke 871 579.56 2986.2 
Strang 721 580.62 2977.5 
Hexagonal 763 573.76 3035.0 
Gmsh 943 579.47 3004.3 
Exact - 579.814 2952 


Ir 


Table 7.4: Table of steady state values for the boundary and divide thickness in the EISMINT 


problem under four different initial meshes. 


number of mesh points increase, the absolute difference between the approximation to 


the boundary and the exact location decreases exponentially. However the amount of 


computational time taken to achieve these results exhibits a linear increase which sug- 


gests there is an optimal trade off between the number of mesh points and the accuracy 


required. We can therefore say the the two-dimensional numerical CMF method is in 


this sense convergent. 
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Figure 7.9: Convergence analysis for varying mesh nodes using the spoked initial mesh: a) 


absolute error in the boundary position, b) computational time. 
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7.5 Summary 


In this chapter we began by looking at the shallow ice PDE in two horizontal spatial 
dimensions. We examined the evolution of the total mass of the domain and described 
the PDE in non-dimensional form with the requirements to maintain balance between 
the terms. 

We then derived the weak formulation for the CMF approach in two dimensions 
which resulted in a net velocity described by the flow velocity and the ice accumulation. 
The weak formulation was then approximated numerically using a finite element method 
to gain the net velocity and a solution. 

The numerical model was then tested against an example scenario that was circular 
in nature, with a zero Dirichlet condition for the ice thickness on the entire boundary. 
This demonstrated that the mesh is able to both advance and retreat whilst maintaining 
connectivity and avoiding mesh tangling. We provided a direct comparison between the 
moving mesh method and the fixed grid approach described in the EISMINT scenario, 
where the results showed a significant increase in accuracy but also in computational 
time. 


Finally we showed the numerical method was convergent. 


Chapter 8 


Data Assimilation on Moving 


Meshes in 2D 


In this chapter we apply the optimised 3D-VAR data assimilation scheme to the two- 
dimensional moving mesh model of the previous chapter to produce a statistically best 
estimate of the state of the ice. We describe the components of the scheme for this 
higher dimensional moving mesh application and highlight the differences with the one- 
dimensional case. We again use a twin experiment to test this approach and demon- 
strate the ability of the method to improve the prediction. 

We then extend the scheme to incorporate the mesh into the state vector as we did 


in the one-dimensional case, highlighting potential restricting factors. 


8.1 3D-VAR on a 2D Moving Domain 


Increasing the numerical model to two horizontal dimensions produces a more accurate 
approximation to the ice thickness over the domain as we allow the ice to spread out 
in more directions. We now apply the 3D-VAR scheme to this moving mesh model 
to further improve the representation of the ice thickness. As in one-dimension the 
2D numerical mesh is also evolving as part of the dynamical model, but this does not 


change the general form of the optimal 3D-VAR equation used before (Eq. (6.1)), given 
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by: 
Zi, = zi + Ki(yen — C,2/), (8.1) 


subject to the dynamical model 

24 = Mazi, (8.2) 
where z; € RR? is a vector containing the unknown state variables at time t,;, where ® 
is the optimal analysis solution, and / is the forecast from the dynamical model. The 


gain matrix K, is given by 
K, = BpC, (C;,B,C; + R,)7!. (8.3) 


There are however, some changes to the observation operator Cz, the observation er- 
ror covariance matrix R,; and the background error covaraince matrix B, in the 2D 


application which we shall detail now. 


8.1.1 Observation Operator 


The observation operator maps the state vector to observational space so that it is 
directly comparable to the observations. 

In one-dimension we defined the rows of the observation operator C; by using a 
simple linear interpolation technique to define an observation of ice thickness in terms 
of the two surrounding mesh points. In two-dimensions this is no longer possible. 
Instead we use the knowledge of the numerical solution which lies on a triangular mesh, 
for which we seek to interpolate an observation to the vertices of a surrounding triangle. 

Consider a point P located at (*,y*) sitting within a triangle with vertices ABC 
(see Figure 8.1). This point may be described in terms of the areal co-ordinates [26] 


La,Lp,Lc, for which 


Apsc 
bA= 8.4 
AE Kise a 
Arca 
ie = 8.5 
3 Aase oe) 
gaa (8.6) 


= ’ 
Asc 
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A 


Figure 8.1: Point P lying within triangle ABC can be expressed in terms of its areal coor- 


dinates. 


where A;;;, represents the area of the triangle with vertices ijk. Since 
La+Lp+Lo=1 (8.7) 


these three co-ordinates provide a simple method of interpolation comparable with the 
linear interpolation used in one-dimension. 
We now write the expressions for [4,02 and Lc in terms of the observation point 
P = (a*,y*). First note that the area of triangle ABC’ with vertices (11,41), (£2, ye), 
(x3, y3) is given by 
dL: ty Gi 
Aagc = ; 1 x2 y2 |: (8.8) 
1 23 Ys 
Equivalent expressions can be formulated for the three smaller triangles in Figure 8.1. 


Therefore L.4 can be written in the form 


1 
La= —(b * 4 day* 8.9 
A 5g iba + can’ + Ay’), (8.9) 
with 
x 1 1 2 
b= a ue > CA= - ; A= : (8.10) 
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Similar expressions for Lg and Lc can be formed with cyclic expressions of bg,cgp,dp,bc,cc 


and dg to provide the interpolations values. The observation may then be written 
A(x*,y") = Lah(x1, y1) + Leh(x2, yo) + Loh(2s, ys) (8.11) 


Therefore for every observation of ice thickness, the corresponding row 7 in C contains 
the coefficients L4,Lg and Lc located at Cj, where 7 represents the numbered location 
of the vertices A,B and C' within the mesh structure. We assume that every observation 


lies within the numerical domain of the model. 


8.1.2 Observation Error Covariance Matrix 


The observation error covariance matrix R, represents the uncertainty in the errors in 
observations and is unaffected by the dimension of the system. As such this matrix can 
be written in the form 


R,=071, R, €R?*?, (8.12) 


under the assumption that each observation is independent of the others. The error 


variance for each observation is given by o?. 


8.1.3. Background Error Covariance Matrix 


The background error covariance matrix B, is a qx q matrix representing the uncertainty 
in the errors in the prior state of the system. 

The elements of the background error covariance matrix B; can be approximated 
using an analytic correlation function, which we choose to be a modified Gaussian 
function. This uses the distance between mesh points to form the elements of B;. For 


two-dimensions (and higher) this may be written as 
bi; = o} exp Elba? Pi ees | (8.13) 


where L is the inverse of a background correlation length scale and o} is the error 


variance associated with each ice thickness variable. The structure of this matrix de- 
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pends on the numbering system used when defining the numerical mesh and like the 
one-dimensional version, the function evolves to reflect this. The background error co- 
variance matrix influences on the way information from observations is spread between 
the variables and by writing it in this form it implicitly depends upon the flow of the 


dynamical model. 


8.1.4 Analysis Solution and Mass Conservation 


Perfoming a step of an assimilation scheme alters the mass of the glacier using external 
information from observations of ice thickness. The CMF moving mesh method that 
we are applying this assimilation scheme to (see Section 7.2.1) is based on conserv- 
ing relative mass, where the conservation principle (Eq. (7.17)) no longer holds after 


assimilating. The conservation principle 


a 
— hdx = pL, 8.14 
A(t) [. ie ory 


used to define the mesh deformation velocity and the ice thickness profile requires an 
analysis solution before the numerical model can continue. 

Therefore at every assimilation cycle, after calculating the analysis solution to the 
ice thickness using Eq. (6.1), we recalculate the total mass 6 such that the analysis 


solution of the mass is 


0“(t) = | h*dx, (8.15) 
Q(t) 
and the mass fraction constants jz are redefined as 
one (On dx = pe (8.16) 
O°(t) Jaw) 


before the forecast step in the dynamical model. The analysis mass and equivalent mass 
fraction constants then correspond to the analysis estimation of the state variables. 
In this work we calculate the analysis solution to the total mass and mass fraction 


constants at the end of each assimilation step. 
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8.1.5 3D-Var Algorithm for the CMF Moving Mesh Method 


We now present the 3D-VAR algorithm for assimilating on a 2D moving domain using 


the weak form conservation of mass fractions method from Section 7.2.1. 


1. Calculate a forecast of the state vector z/ by evolving the numerical model of the 
CMF method (Section 7.3) using the previous analysis solution as model inputs. 
This gives forecast values of the numerical mesh (X, Y)/, the ice thickness H/ 


and the total mass 6/. We also know the constant in time mass fractions ju. 


2. Use the optimal 3D-VAR formulation (Eq. (8.1)) to produce a new analysis solu- 
tion of the ice thickness H® using observations available at the current time. The 


analysis solution to the numerical mesh is identical to the forecast, with X* = X/. 


3. Calculate the analysis solution of the total mass 6% and update the mass fractions 


* using Eqs. (8.15) and (8.16). 


4. Evolve the analysis solution using the dynamical model to produce a forecast at 


the next time step. 


5. Repeat steps 2-4. 


8.2 A Test Experiment 


We now test the 2D data assimilation scheme with the components defined above using 
a test scenario. As there are no observations available we will use a twin experiment, 
as in the 1D method (Section 6.2). A reference solution is simulated by evolving the 
numerical model in Section 7.4.1 with exact initial conditions, from which we sample 
observations. We then define incorrect initial conditions for the forecast model which 


we will evolve and apply the data assimilation scheme to. 
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8.2.1 Initial Conditions 
Reference Solution 


To generate the observations for a twin experiment we first define the known set of 


exact initial conditions. Let the true initial ice thickness be the function 
h=1-(¢+y’), 2 +yY <1. (8.17) 


The model is evolved using the finite element approximation described in Section 7.3, 
with At = 0.2a on a Hexagonal-type mesh (see Appendix A.3) with 187 mesh nodes. 
Both the reference and forecast models run until an end time T’ = 2000a, as we want to 
see the effects of including observations before the models reach a steady state solution. 
These initial conditions formulate the reference solution, which provide observations 


and the target solution. 


Forecast Model 


For our experiments we estimate the initial ice thickness by function 
h=12—-(2?4+y*), 2 t+y? <1.2. (8.18) 


where the estimate of the initial domain, in addition to the ice thickness profile, is 
different to the reference solution. The total mass will also be different, along with the 
mass fraction constants used in the moving mesh procedure. The time step is taken to 
be the same as the reference solution, and we again use a Hexagonal-type mesh with 


187 mesh nodes. 


8.2.2 Results 


We include a set of 18 observations of ice thickness taken at t = 500 and t = 1500, 
distributed randomly across the domain with random noise taken from a normal dis- 


tribution ¢, ~ N(0,02) with o? = 0.05 to represent an error of 5%. The observation 
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noise is kept constant between each of the experiments. The background error variance 
in Eq. (8.13), o? is also set to 0.05. 

We see in Figures 8.2 and 8.3 that the application of the 3D-VAR scheme to the 
two-dimensional moving mesh method provides improved results for an ice thickness 
profile, with an small improvement to the extent of the domain. 

At the first assimilation time, t = 500, that there is a large difference between the 
forecast (Figure 8.2(a)) and reference (Figure 8.2(c)) thickness profiles. After the assim- 
ilation the analysis thickness profile (Figure 8.2(b)) more closely resembles the truth, 
particularly in the centre of the domain. When we examine the difference between 
the analysis and reference solutions, Figure 8.2(d), there is a notable ring relating to 
the largest difference that corresponds to the location of the boundary of the reference 
solution. This is because the shape and size of the analysis domain remains the same 
size as the forecast since we only update the ice thickness variables not the domain. 

By the second assimilation time, t = 1500, we see the see the forecast domain 
(Figure 8.3(a)) is closer to the reference (Figure 8.3(c)), as the moving mesh method 
implicitly improves the estimate of the mesh as it evolves. A second assimilation step 
again improves the prediction of the ice thickness profile (Figure 8.3(b)), with a smaller 
difference between the reference and the analysis as seen in Figure 8.3(d) which again 
exhibits the largest difference where the reference boundary is located. 

The evolution of one of the boundary nodes is shown in Figure 8.4 where we can 
directly observe this moving point being pulled towards the reference location after 
the first assimilation time despite no analysis correction to the boundary at that time. 
Interestingly, after the second assimilation time the boundary appears to be moving 
with the same speed as the reference boundary. This suggests that the ice thickness 
profile is similar enough to the reference solution to generate the correct mesh velocity, 
despite the domain being overestimated. 

Similar results are achieved using the alternative initial meshes described in Ap- 
pendix A, which have little effect on the overall impact of the observations shown 
already. Additionally, experiments exploring the impact of observation quantity and 


distribution consistently improve the overall accuracy of the scheme. For instance, in- 
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Forecast Analysis 


Reference 


Figure 8.2: 2D moving mesh with data assimilation at t = 500:- a) Forecast before assim- 
ilation, b) Analysis solution, c) Reference solution and d) Difference between analysis and 


reference. 


cluding observations all on one side of the domain instead of distributed randomly does 
not significantly impact the symmetry in the dynamical model. 

Observations that lie within elements connected to the boundary nodes occasionally 
caused spurious results, which are believed to be due to the steep gradients along with 


the enforcement of Dirichlet boundary conditions. 
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Figure 8.3: 


similation, b) Analysis solution, 
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2D moving mesh with data assimilation at t = 1500:- a) Forecast before as- 


c) Reference solution and d) Difference between analysis and 
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2D moving mesh with data assimilation:- Evolution of a boundary point over 


time- Forecast (blue), Reference (Green), Analysis (Red). 
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8.3 2D Extended Scheme: Assimilating the Mesh - 


Theory 


In Section 6.3 we built an extension to the one-dimensional 3D-Var scheme to incor- 
porate the mesh explicitly into the state vector to calculate an analysis solution to the 
numerical domain. This was successfully implemented and it was demonstrated that 
the prediction of the domain at the time of assimilation was improved. We now seek a 
similar extension in higher spatial dimensions. 

Generally, positional features in higher dimensions are observed in the form of a 
line, such as the boundary of the domain or geological evidence such as moraines for 
glaciers. The interpolation of this to a numerical mesh requires a discretisation of the 
line into points that may be related to the mesh via the operator Cy. 

The definition of this operator, and the error covariance matrices R; and By is not 
as straightforward in higher dimensions, since each co-ordinate is now a vector. We 


present two options for solving the problem in 2D: 


8.3.1 Option 1: Conjoined Co-ordinates 


The first option joins the two co-ordinates together to form a single entity, d; = ||x;,||, 


such that the state vector is 


i= (8.19) 
h 


By combining the co-ordinates in this way, we essentially reduce the problem to 1D 
and can define the other components in the data assimilation scheme in a similar way 
to Section 6.3. The extended observation operator to map the state vector z. to the 


observations can be written as 


Cu 0 
C.= ; (8.20) 
0 C 


where Cy, is a simple interpolation matrix mapping observations of positional locations 


(e.g. the moving boundary) to the equivalent numerical mesh points. 
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The extended observation error covariance matrix, assuming observations are inde- 


pendent is 


o2,1 O 
R. = (8.21) 
0 R 


where o2, represents the error variances relating to observations of positional features. 


Lastly, the background error covariance matrix can be written in the form 


Bu B 
B. = ee ails (8.22) 
(Bu)? B 


where Bag is the error covariance matrix corresponding to the numerical mesh. As in 
one-dimension this can be defined using the modified Gaussian function from Eq. (8.13) 
as a correlation function. The cross-covariances Bg, may then also be determined using 


the boundary points as was done successfully in one-dimension. This approach means 


the matrix B, € R?%*?% is straightforward to calculate and suffers from no increase in 


size over the one-dimensional method. 
With this option the analysis step, Eq. (8.1), is no more difficult than the one-dimensional 
method. However, by treating the coordinate pair as a single item, some post-assimilation 


technique is required to distinguish between the x and y elements. 


8.3.2 Option 2: Separate Co-ordinates 


An alternative option separates the co-ordinate pairs and treats them as separate vari- 


ables, such that the state vector is defined by 


with x and y representing the sets of x and y co-ordinates respectively. The components 
of the data assimilation algorithm can then be determined as follows. 
By separating the variables any observations of positional features, which contain 


information in the form of a co-ordinate pair (xz*,y*), will also be separated. The 
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extended observation operator to map the state vector to these observations will be of 


the form: 

C, O O 
Cy 0 
0 C 


C= (8.24) 


0 
0 
where C,, interpolates the x co-ordinates in the state vector to the observation location 
x* and C, interpolates the y co-ordinates to y*. These matrices need to account for 
the net distance between the observation and the mesh points rather than treating the 
co-ordinates completely separate since, for example, locations that are close together 
in the x-direction could be far apart in the y-direction. 

Assuming that errors in observations are independent of each other, the extended 


observation error covariance matrix may be written in the form 


0 oO R 
Re=]|o7,1 0 0 (8.25) 
0 «2,1 0 


where o2, and @;,, represent the error variances relating to observations of positional 
features in the x and y cartesian directions. We assume that the error variances are the 
same, 02, = 93,. The assumption that errors in the observations of positional features 
are independent is debatable since we are observing a co-ordinate pair. 
Finally, the extended background error covariance matrix for the error covariances 
in the forecast can be written as 
Boz Bry Ban 
Be = (B,,)" By — Byn (8.26) 
(Ben)’ (By)" B 
where B,, and B,, are the error covariance matrices for the forecast of the x and y 
co-ordinates, while B,,, Bz, and B,, are the cross-covariance matrices between the 
co-ordinates and the ice thickness. 
Defining Braz, Biz, Ban,Bry,Brn and By» is certainly a grey area. As previously 
mentioned the co-ordinates cannot be treated completely separately as the net distance 


between mesh points needs to be considered. 
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8.3.3 Discussion 


Option 1 joined the co-ordinates together into a single unit, which reduces the problem 
to one-dimension and makes it easier to define and calculate. In doing so we require a 
method to distribute the analysis solution back to the individual x and y components. 

Option 2 splits the co-ordinate pairs into separate entities, which maintains the 
full description of all the variables but is more difficult to define. The construction 
of the components in the data assimilation scheme still needs consider the distance 
between mesh points when calculating the observation operator mapping and the error 
covariances. 

In essence, both methods appear to have significant drawbacks but until they have 
been tested we cannot determine which approach, if either, provides an viable analysis 


solution to the numerical mesh. 


8.4 Summary 


In this chapter we presented an application of the 3D-VAR data assimilation scheme 
to a two-dimensional moving mesh environment. An analysis solution of the total and 
relative mass values was required to ensure compatibility with the chosen method of 
mesh movement. Results demonstrated a successful application of the scheme with 
improved predictions of the ice thickness profile. 

We then attempted to extend the scheme to include the mesh within the state vector 
and allow the potential for observations of positional co-ordinates to be included, using 
a similar approach to the extension in one-dimension. However, since the mesh nodes 
are described in coordinate pairs this proved significantly more difficult. A couple of 
potential solutions was proposed, neither of which is ideal. This remains an open ended 


problem. 


Chapter 9 


Conclusions 


Accurate, efficient numerical modelling of the dynamical flow of glaciers is vital to 
simulating and forecasting the Cryosphere, which in turn has an impact on the atmo- 
sphere and oceans of our planet. Using adaptive mesh techniques improve efficiency 
by reducing the need for numerical domains with high resolutions and in particular 
moving meshes adapt the mesh by relocating mesh points to areas that require higher 
resolution. These numerical models can make use of the global observation network to 
produce a statistically best representation of glaciers using data assimilation. 

In this thesis we sought to apply both a moving mesh method and data assimilation 
to the shallow ice model used to simulate glacier flow. We now summarise the work 
undertaken in this thesis before drawing conclusions and discussing ideas for future 


research. 


9.1 Summary 
In the introduction we set out the objectives of this thesis. They were: 


e Develop a moving mesh method to equations governing dynamical ice flow. 


e Use data assimilation to combine observations with a moving mesh model of 


dynamical ice flow. 
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With the specific aims of: 


e Applying the CMF moving mesh method of Baines, Hubbard and Jimack [8] to 


the one-dimensional shallow ice approximation equation used within glaciology. 
e Analysing the impact of internal glacier flow on the movement of the ice boundary. 


e Using a sequential 3D-Var data assimilation scheme in conjuction with the moving 


mesh to improve the predicted ice profile. 


e Building an extended scheme capable of including observations of positional fea- 
tures within the data assimilation method to improve the representation of the 


domain. 


e Extending the moving mesh method and the 3D-Var scheme to two horizontal 
dimensions and demonstrate the difficulties encountered when considering this 


more complex scenario. 


In the introduction we discussed elements of Glaciology and the reasons for being 
interested in ice dynamics. We continued this in Chapter 2 by introducing some of 
the terminology used in the field and the key driving forces behind ice flow. A brief 
discussion of the types of models used to simulate this flow was given, with particular 


focus on the shallow ice equations. 


Chapter 3 introduced the idea of numerical models as a solution to time depen- 
dent PDEs, focusing on adaptive mesh models. We discussed the different methods of 
adaptivity, with the benefits of each method for ice flow models presented. A detailed 
description of a particular moving mesh method [8] was provided, based upon the idea 
of conserving relative mass. We then provided an brief history of the use of adaptive 


mesh methods within the field of glaciology. 


In chapter 4 we introduced the idea and objectives of data assimilation, with an 


overview of some of these methods. Both sequential and variational schemes were 
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introduced, with details of the various components that are required to apply data 
assimilation techniques. We then examined the past attempts at applying data assim- 


ilation to both dynamical ice sheet models and adaptive mesh methods. 


We began chapter 5 by examining the one-dimensional shallow ice equation in more 
detail. In particular we evaluated the effect of internal diffusive flow on the movement 
at the boundary. This provided a criterion upon the local ice profile that allows the 
boundary to move in the absence of any external forcing. We introduced a method for 
calculating a steady state solution when the forcing term is independent of time. We 
proceeded to apply the relative mass conservation method to the shallow ice equation, 
relating the resulting expression for mesh movement back to the analysis of boundary 
movement. A finite difference approximation to this moving mesh method was devel- 
oped and tested using a series of simple test equations. We also extended the method 
to allow for radially symmetric problems and compared the method to the circular EIS- 


MINT moving margin experiment. 


Chapter 6 detailed the application of the 3D-Var data assimilation scheme to the 
1D moving mesh method of the previous chapter. We examined the changes to the 
components that make up the data assimilation scheme for a moving mesh method and 
described the analysis update of the total mass and the relative mass constants that 
were required in order to apply the assimilation scheme to our method. We then tested 
this method using a twin experiment applied to the test equations introduced in the 
previous chapter. We made an extension to the data assimilation scheme to include 
the mesh within the state vector and allow observations of positional features to be 
included. The impact on the components was assessed, in particular the background 
error covariance matrix with a proposed definition of the cross covariances between 
the mesh and the ice thickness variables. We carried out experiments to demonstrate 
the impact of including observations of the boundary and cross covariance terms, both 


separately and simultaneously. 
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In chapter 7 we extended the shallow ice equation to two horizontal dimensions, 
and detailed a non-dimensional version of the equations. The moving mesh method 
was then described in a weak formulation, from which we extracted the expression for 
the mesh velocity. A finite element approximation to the moving mesh method was ap- 
plied and tested using a set of test scenarios similar to those used in one-dimension. We 
then compared the method to the EISMINT moving margin experiment and assessed 


the impact of using different initial meshes. 


Lastly, in chapter 8 we applied the 3D-Var scheme to the 2D moving mesh model 
and detailed changes to the components of the scheme when an additional horizontal 
dimension is included. We used a twin experiment to test the scheme against the test 
equations and the EISMINT scenario. Finally we attempted to apply an extension 
to the scheme to include the mesh in the state vector and detailed the difficulties 


encountered. 


9.2 Conclusions 


In conclusion, we have shown that moving mesh methods accurately simulate the dy- 
namical flow of glaciers and can achieve greater accuracy than conventional fixed grid 
approaches for, in 1D at least, minimal computational cost. A 3D-VAR Data assimila- 
tion method can be applied to moving mesh methods after adjustmenting the compo- 
nents of the scheme and can, in 1D at least, be extended to incorporate the numerical 
mesh within the analysis update. 


Specifically, we conclude that 


e The CMF moving mesh method of Baines, Hubbard and Jimack [8] is an effective 
tool for simulating the dynamical flow of a one-dimensional shallow ice flowline 


model. 
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We demonstrated in chapter 5 that the moving mesh approach can effectively 
simulate various aspects of dynamical ice flow in one-dimension. This method uses 
a combination of the diffusive velocity and the ice equivalent accumulation rate to 
adapt the mesh and accurately simulate changes in the domain. Results show that 
when using the same number of grid points, the solution to the EISMINT moving 
margin experiment, especially the location of the boundary, is significantly closer 


to the true solution for a negligible computational cost. 


e Internal glacier flow impacts the movement of the ice boundary when the local 


profile has an asymptotically infinite gradient. 


In Section 5.2.1 we assessed the diffusive flow velocity by rewriting the equation 
in a different form. This led to a condition on the local ice thickness profile that 
needed to be met to induce a finite flow at the boundary, corresponding to an 
infinite gradient. The underlying topography, ice density and temperature had no 
impact on this condition. This was supported by experiments with a zero source 


term at the initial boundary. 


e A sequential 3D-Var data assimilation scheme can be used in conjuction with the 


moving mesh to improve the accuracy of the predicted ice profile. 


In chapter 6 we demonstrated that with little configuration a 3D-Var scheme is 
readily applicable to the moving mesh scenario. An analysis update of the total 
mass and the relative mass constants was required for this particular moving 
mesh scheme to be compatible with the data assimilation. A modified Gaussian 
function was used to form the background error covariance matrix, which gave 
an improved representation of the evolving mesh that was previously unchanged 
when assimilating. Results showed that by observing the ice thickness the profile 
is improved, which with future iterations of the numerical model also reduces the 


errors in the domain. 
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e An extended scheme allows an analysis solution to the mesh at the time of assim- 


ilation. 


By explicitly including the mesh into the state vector in Section 6.3, we were able 
to calculate an analysis solution of the numerical mesh. This was achieved in two 
ways; firstly by allowing observations of physical points such as the boundary and 
secondly by defining the cross-covariance terms in the background error covariance 
matrix to pass information between the mesh and the ice thickness variables. The 
experiments show that including one or both of these results in an improvement 
to the prediction of the numerical domain in addition to the ice thickness profile 


at the moment of assimilation. 


e The moving mesh method readily extends to two horizontal dimensions. 


We showed in chapter 7 that rewriting the moving mesh method in a weak for- 
mulation allows us to use finite elements to approximate the solution to the shal- 
low ice equation in two-dimensions. Experiments demonstrate the ability of the 
method to replicate the different types of movement. Comparisons with the EIS- 
MINT moving margin scenario show that using erds the number of mesh nodes 
as a fixed grid method achieves a significant reduction in error, particularly the 
approximation to the location of the boundary. The choice of initial mesh demon- 


strated small variability in the steady state solution achieved. 


e The sequential 3D-Var scheme applied to the moving mesh method extends to 


two horizontal dimensions. 


In chapter 8 we demonstrated that the 3D-Var scheme can be applied to the 
moving mesh method in two-dimensions with a few modifications to the schemes 
components. An analysis update of the total mass and the relative mass con- 


stants was required for compatibility. Results indicated an improvement in the 
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prediction of the ice profile, with the error in the domain representation reducing 


implicitly with iterations of the numerical model. 


e Extending the 3D-Var scheme to achieve an analysis solution of the mesh presents 


many difficulties in two-dimensions. 


Two options for explicitly including the numerical mesh in the data assimilation 
scheme in 2D where proposed, with their advantages and disadvantages discussed. 


The effectiveness of both options is yet to be tested. 


9.3. Future Work 


In this work we have used the shallow ice approximation which defines an explicit equa- 
tion to represent the diffusive flow velocity, where the resistance is solely at the bed 
of the glacier. In other scenarios different resistive stresses become important. For 
instance, in ice streams the resistance at the edge of the moving ice is more important 
than the base, while in problems that involve an ice shelf the back stress is most critical. 
In theory any time-independent function to describe the flow velocity can readily be 
included within the moving mesh solution to the mass continuity PDE (Section 5.3.1), 
which means the moving mesh procedure is not restricted to the shallow ice approxi- 
mation. However, solving the function for diffusive flow is a separate issue independent 


of the moving mesh. 


In our experiments we neglected the influence of temperature and density by includ- 
ing them as constant parameters in the model. Temperature can be incorporated by 
solving a thermodynamic equation, which in turn influences the flow law parameter A 
in Eq. (2.3). Including variable density is a different issue, as the main changes occur 
with depth. Potentially, a multilayered model could be a solution which would also 
allow an element of flow in vertical directions; however, the form of the mass continuity 


equation would change. A multilayered moving mesh solution would need careful defi- 
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nition around the boundary layer, but would open up a whole new area of development. 


There are other scenarios in glaciology that are of particular interest, most notably 
the issue of finding the grounding line which separates grounded ice from floating ice 
shelves. With a moving mesh model we could fix a mesh node to the grounding line 
itself which would enable us to track its movements rather than interpolating to find 
its location. Additionally, if we were able to observe the grounding line the extended 


data assimilation scheme (Section 6.3) would be ideal for including this observation. 


The glacier bed has been assumed to be flat throughout this thesis, except in the 
analysis of the diffusive velocity at the boundary (Section 5.2.1). The choice of bed 
will have an effect on the diffusive flow, particlarly in areas where the gradient of the 
bed is positive. In addition more complex models aim to include the basal sliding that 
glaciers exhibit when there is water present at the base. Furthermore, over long time 
scales the effect of isostasy will influence the shape of the bed. The inclusion of the 
glacier bed opens up the possibility to apply the moving mesh method to a real-life 
scenario, which provides the opportunity to fully assess the feasibility of the method to 


accurately simulate dynamical flow. 


In our data assimilation model the biggest limitation involved assuming all observa- 
tions occur within the computational domain. The next logical step would be to devise 
a method for incorporating observations from outside the domain, ideally into the ex- 
tended scheme so that the numerical mesh is corrected given this new information. The 
next step for the extended scheme is to find an effective way of applying the method to 
higher dimensions, either by developing the ideas mentioned in Section 8.3 or devising 


an alternative method. 


When the model has been applied to a real-life scenario, actual observations can 


be included into the data assimilation scheme to further test the abilities of both the 
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original and extended schemes to improve the estimate of the ice thickness profile and 
numerical domain. Including observations of the diffusive flow (at the surface) opens 
up the possibility of rewriting the problem by expressing the state vector in terms of 
this flow instead of the ice thickness. There could then be two applications of a data 
assimilation scheme at each time step; once to gain an analysis estimation of the dif- 
fusive flow before the mesh is moved and once after to update the ice thickness (and 


domain if positional observations are available). 


Long term, applying an ensemble data assimation scheme (see e.g. [32]) to the mov- 
ing mesh environment is a key objective. These methods use an ensemble of perturbed 
initial conditions, forecast with the dynamical model to estimate the errors in the model 
and generate the optimal forecast. In a moving mesh method the ensemble approach 


would provide the best possible estimate of the numerical domain. 


Appendices 
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Appendix A 


Circular Mesh Generation 


There are many ways to discretise a domain. Here we show four methods for discretising 


a circular shaped domain around the origin. 


A.1 Spoke 


One of the most straightforward circular meshes is a spoked mesh, with the nodes lying 
along straight lines extending from the centre. like a bike wheel. Two parameters are 
used to determine the number of grid points; the number of points on each spoke and 
the number of spokes themselves. If we define the number of spokes in terms of the 
angle between them, this gives a set of polar co-ordinates with which to generate the 
mesh. 

Whilst simple to calculate, as you can see in Figure A.1 the nodes are concentrated 
close to the centre and spread out around the boundary. This would be better suited 
to a problem where the centre requires more information, as opposed to the ice sheet 


problem where the boundary requires a finer mesh. 
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Figure A.1: An example of a mesh generated using spokes pertruding from the centre, with 


5 nodes per spoke and the angle between them 7/8. 
A.2 Strang 


An alternative mesh is the one proposed by Strang [98]. Here the nodes are placed 
using a combination of spokes and circles. The nodes still lie on the spokes however for 
each circle there is only a node on every other spoke. The next circle out then places a 
node on remaining spokes and then the process is repeated. The number of circles and 
the number of nodes lying on each circle then define the mesh. 

This mesh is relatively simple to produce and around the centre there is a smaller 
cluster of nodes than the spoke mesh. However there are still large connections between 
nodes near the boundary and long thin triangles are produced from the calculation (see 


Figure A.2). This effect is reduced by including more nodes. 
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Figure A.2: A mesh generated using the Strang approach with 10 circles and 8 nodes per 
circle. At low resolution this mesh loses information close to the boundary and has a lot of 


long thin triangles which can cause numerical problems due to the shallow angle. 
A.3 Hexagonal 


Looking closely at Figure A.1(a) and Figure A.2(a) shows that these two methods of 
mesh generation can lead to nodes with varying amounts of connections. With the 
exception of the boundary, it can be beneficial to have an equal number of connections 
to improve the structure of the matrices in the finite element method Section 7.3. With 
this in mind we construct a mesh where the elements are hexagonal in shape with 
six connections into each node. To achieve this we again define how many circles we 
require, however this time for each circle we double the number of nodes each time as 
we move away from the centre, starting from six. The nodes are then distributed evenly 
around the circle. 

As we observe in Figure A.3 the interior nodes now contain six connections and there is 
a greater number of nodes near the boundary to provide more detail, but this introduces 


long thin triangles. 
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Figure A.3: A hexagonal mesh with only 4 circles. Note that each interior node has 6 


connections each and the regular hexagon in the centre. 


A.4 Automatic Generation 


The final mesh is one that is automatically generated using the open source software 
GMSH [37]. The software discretises along the lines of the shape, in this case the 
circumference and diameter of the circle, which are then refined to produce a mesh. 


Further refinement can be performed to achieve the level of resolution required. 
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Figure A.4: Mesh generated using the software GMSH [87]. The software uses a refinement 


criteria that leaves the individual triangular elements approximately equilateral in shape. 


Appendix B 


Useful Theorems 


We state some of the theorems used in this thesis for convenience, the proofs for which 


can be readily found. 


B.1  Leibniz’s Integral Rule 


The differentiation of a definite integral, with limits that are functions of the differential 


variable is given by 


do fo 


a 
inf omen i £ fle,ude + Holy) 


B.2 L’Hopital’s Rule 


L’Hopital’s rule states that if 


lim f(a) = lim g(x) = O(orce), 


wa wa 


then 
lim f(z) = lim ; 
esa g(r) 2a g(x) 


provided the derivatives exist. 


148 


(B.1) 
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B.3 Trapezium Rule 


An integral may be approximated using the Trapezium Rule, where the region is split 


into N mesh points and approximated by 
[1 fle 


B.4 Reynold’s Transport Theorem 


Sto &i)(f(tisa) + f(#a)). (B.4) 


ae 


A general three dimensional version of Leibniz’s Integral Rule is given by Reynolds 


Transport Theorem [90]: 


d 


d 
a ge t)dx = [. qt & Diet f (wan) t)dT, (B.5) 


where n is the unit normal vector and v is the velocity of the area element. 


B.5 Divergence Theorem 


The Divergence Theorem relates the volume integral of the divergence to the flux on 


the surface in the outward normal direction. This is written as: 


[ev seoyax= $ (r60).myar (B.6) 
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